MAINT: Refactor polynomial leastsquares fits.

* Rewrite the documentation for `deg`.
* Check that the passed degrees are int and raise if not.
* Do not accept `deg` arguments that do not convert as 0-D or 1-D.
* Sort passed degree terms for repeatability.

7 files changed, 105 insertions, 138 deletions

diff --git a/numpy/polynomial/_polybase.py b/numpy/polynomial/_polybase.py
index 37eb59f03..6fa72b6f9 100644
--- a/numpy/polynomial/_polybase.py
+++ b/numpy/polynomial/_polybase.py
@@ -739,14 +739,11 @@ class ABCPolyBase(object):
             y-coordinates of the sample points. Several data sets of sample
             points sharing the same x-coordinates can be fitted at once by
             passing in a 2D-array that contains one dataset per column.
-        deg : int or array_like
-            Degree of the fitting polynomial. If `deg` is a single integer
-            all terms up to and including the `deg`'th term are included.
-            `deg` may alternatively be a list or array specifying which
-            terms in the Legendre expansion to include in the fit.
-
-            .. versionchanged:: 1.11.0
-            `deg` may be a list specifying which terms to fit
+        deg : int or 1-D array_like
+            Degree(s) of the fitting polynomials. If `deg` is a single integer
+            all terms up to and including the `deg`'th term are included in the
+            fit. For Numpy versions >= 1.11 a list of integers specifying the
+            degrees of the terms to include may be used instead.
         domain : {None, [beg, end], []}, optional
             Domain to use for the returned series. If ``None``,
             then a minimal domain that covers the points `x` is chosen.  If
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 9db613b78..2537bea32 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -1617,14 +1617,11 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int or array_like
-        Degree of the fitting series. If `deg` is a single integer
-        all terms up to and including the `deg`'th term are included.
-        `deg` may alternatively be a list or array specifying which
-        terms in the Legendre expansion to include in the fit.
-
-        .. versionchanged:: 1.11.0
-        `deg` may be a list specifying which terms to fit
+    deg : int or 1-D array_like
+        Degree(s) of the fitting polynomials. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included in the
+        fit. For Numpy versions >= 1.11 a list of integers specifying the
+        degrees of the terms to include may be used instead.
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1718,11 +1715,11 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     """
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
-    deg = np.asarray([deg,], dtype=int).flatten()
+    deg = np.asarray(deg)
 
     # check arguments.
-    if deg.size < 1:
-        raise TypeError("expected deg to be one or more integers")
+    if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+        raise TypeError("deg must be an int or non-empty 1-D array of int")
     if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
@@ -1734,19 +1731,17 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
-    if deg.size == 1:
-        restricted_fit = False
-        lmax = deg[0]
+    if deg.ndim == 0:
+        lmax = deg
         order = lmax + 1
+        van = chebvander(x, lmax)
     else:
-        restricted_fit = True
-        lmax = deg.max()
-        order = deg.size
+        deg = np.sort(deg)
+        lmax = deg[-1]
+        order = len(deg)
+        van = chebvander(x, lmax)[:, deg]
 
     # set up the least squares matrices in transposed form
-    van = chebvander(x, lmax)
-    if restricted_fit:
-        van = van[:, deg]
     lhs = van.T
     rhs = y.T
     if w is not None:
@@ -1776,11 +1771,11 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     c = (c.T/scl).T
 
     # Expand c to include non-fitted coefficients which are set to zero
-    if restricted_fit:
+    if deg.ndim > 0:
         if c.ndim == 2:
-            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+            cc = np.zeros((lmax + 1, c.shape[1]), dtype=c.dtype)
         else:
-            cc = np.zeros(lmax+1, dtype=c.dtype)
+            cc = np.zeros(lmax + 1, dtype=c.dtype)
         cc[deg] = c
         c = cc
 
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 5d4b357fe..e234c8e23 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -1388,14 +1388,11 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int or array_like
-        Degree of the fitting polynomial. If `deg` is a single integer
-        all terms up to and including the `deg`'th term are included.
-        `deg` may alternatively be a list or array specifying which
-        terms in the Legendre expansion to include in the fit.
-
-        .. versionchanged:: 1.11.0
-        `deg` may be a list specifying which terms to fit
+    deg : int or 1-D array_like
+        Degree(s) of the fitting polynomials. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included in the
+        fit. For Numpy versions >= 1.11 a list of integers specifying the
+        degrees of the terms to include may be used instead.
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1494,11 +1491,11 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     """
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
-    deg = np.asarray([deg,], dtype=int).flatten()
+    deg = np.asarray(deg)
 
     # check arguments.
-    if deg.size < 1:
-        raise TypeError("expected deg to be one or more integers")
+    if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+        raise TypeError("deg must be an int or non-empty 1-D array of int")
     if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
@@ -1510,19 +1507,17 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
-    if deg.size == 1:
-        restricted_fit = False
-        lmax = deg[0]
+    if deg.ndim == 0:
+        lmax = deg
         order = lmax + 1
+        van = hermvander(x, lmax)
     else:
-        restricted_fit = True
-        lmax = deg.max()
-        order = deg.size
+        deg = np.sort(deg)
+        lmax = deg[-1]
+        order = len(deg)
+        van = hermvander(x, lmax)[:, deg]
 
     # set up the least squares matrices in transposed form
-    van = hermvander(x, lmax)
-    if restricted_fit:
-        van = van[:, deg]
     lhs = van.T
     rhs = y.T
     if w is not None:
@@ -1552,7 +1547,7 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     c = (c.T/scl).T
 
     # Expand c to include non-fitted coefficients which are set to zero
-    if restricted_fit:
+    if deg.ndim > 0:
         if c.ndim == 2:
             cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
         else:
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index da441af83..08e83899a 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -1385,14 +1385,11 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int or array_like
-        Degree of the fitting polynomial. If `deg` is a single integer
-        all terms up to and including the `deg`'th term are included.
-        `deg` may alternatively be a list or array specifying which
-        terms in the Legendre expansion to include in the fit.
-
-        .. versionchanged:: 1.11.0
-        `deg` may be a list specifying which terms to fit
+    deg : int or 1-D array_like
+        Degree(s) of the fitting polynomials. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included in the
+        fit. For Numpy versions >= 1.11 a list of integers specifying the
+        degrees of the terms to include may be used instead.
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1491,11 +1488,11 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     """
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
-    deg = np.asarray([deg,], dtype=int).flatten()
+    deg = np.asarray(deg)
 
     # check arguments.
-    if deg.size < 1:
-        raise TypeError("expected deg to be one or more integers")
+    if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+        raise TypeError("deg must be an int or non-empty 1-D array of int")
     if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
@@ -1507,19 +1504,17 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
-    if deg.size == 1:
-        restricted_fit = False
-        lmax = deg[0]
+    if deg.ndim == 0:
+        lmax = deg
         order = lmax + 1
+        van = hermevander(x, lmax)
     else:
-        restricted_fit = True
-        lmax = deg.max()
-        order = deg.size
+        deg = np.sort(deg)
+        lmax = deg[-1]
+        order = len(deg)
+        van = hermevander(x, lmax)[:, deg]
 
     # set up the least squares matrices in transposed form
-    van = hermevander(x, lmax)
-    if restricted_fit:
-        van = van[:, deg]
     lhs = van.T
     rhs = y.T
     if w is not None:
@@ -1549,7 +1544,7 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     c = (c.T/scl).T
 
     # Expand c to include non-fitted coefficients which are set to zero
-    if restricted_fit:
+    if deg.ndim > 0:
         if c.ndim == 2:
             cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
         else:
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 280e28159..d459551ae 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -1387,14 +1387,11 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int or array_like
-        Degree of the fitting polynomial. If `deg` is a single integer
-        all terms up to and including the `deg`'th term are included.
-        `deg` may alternatively be a list or array specifying which
-        terms in the Legendre expansion to include in the fit.
-
-        .. versionchanged:: 1.11.0
-        `deg` may be a list specifying which terms to fit
+    deg : int or 1-D array_like
+        Degree(s) of the fitting polynomials. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included in the
+        fit. For Numpy versions >= 1.11 a list of integers specifying the
+        degrees of the terms to include may be used instead.
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1493,11 +1490,11 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     """
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
-    deg = np.asarray([deg,], dtype=int).flatten()
+    deg = np.asarray(deg)
 
     # check arguments.
-    if deg.size < 1:
-        raise TypeError("expected deg to be one or more integers")
+    if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+        raise TypeError("deg must be an int or non-empty 1-D array of int")
     if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
@@ -1509,19 +1506,17 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
-    if deg.size == 1:
-        restricted_fit = False
-        lmax = deg[0]
+    if deg.ndim == 0:
+        lmax = deg
         order = lmax + 1
+        van = lagvander(x, lmax)
     else:
-        restricted_fit = True
-        lmax = deg.max()
-        order = deg.size
+        deg = np.sort(deg)
+        lmax = deg[-1]
+        order = len(deg)
+        van = lagvander(x, lmax)[:, deg]
 
     # set up the least squares matrices in transposed form
-    van = lagvander(x, lmax)
-    if restricted_fit:
-        van = van[:, deg]
     lhs = van.T
     rhs = y.T
     if w is not None:
@@ -1551,7 +1546,7 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     c = (c.T/scl).T
 
     # Expand c to include non-fitted coefficients which are set to zero
-    if restricted_fit:
+    if deg.ndim > 0:
         if c.ndim == 2:
             cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
         else:
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 2035ba6e9..54e9895db 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -1418,14 +1418,11 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int or array_like
-        Degree of the fitting polynomial. If `deg` is a single integer
-        all terms up to and including the `deg`'th term are included.
-        `deg` may alternatively be a list or array specifying which
-        terms in the Legendre expansion to include in the fit.
-
-        .. versionchanged:: 1.11.0
-        `deg` may be a list specifying which terms to fit
+    deg : int or 1-D array_like
+        Degree(s) of the fitting polynomials. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included in the
+        fit. For Numpy versions >= 1.11 a list of integers specifying the
+        degrees of the terms to include may be used instead.
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1521,11 +1518,11 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     """
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
-    deg = np.asarray([deg,], dtype=int).flatten()
+    deg = np.asarray(deg)
 
     # check arguments.
-    if deg.size < 1:
-        raise TypeError("expected deg to be one or more integers")
+    if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+        raise TypeError("deg must be an int or non-empty 1-D array of int")
     if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
@@ -1537,19 +1534,17 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
-    if deg.size == 1:
-        restricted_fit = False
-        lmax = deg[0]
+    if deg.ndim == 0:
+        lmax = deg
         order = lmax + 1
+        van = legvander(x, lmax)
     else:
-        restricted_fit = True
-        lmax = deg.max()
-        order = deg.size
+        deg = np.sort(deg)
+        lmax = deg[-1]
+        order = len(deg)
+        van = legvander(x, lmax)[:, deg]
 
     # set up the least squares matrices in transposed form
-    van = legvander(x, lmax)
-    if restricted_fit:
-        van = van[:, deg]
     lhs = van.T
     rhs = y.T
     if w is not None:
@@ -1579,7 +1574,7 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     c = (c.T/scl).T
 
     # Expand c to include non-fitted coefficients which are set to zero
-    if restricted_fit:
+    if deg.ndim > 0:
         if c.ndim == 2:
             cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
         else:
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 7c922c11b..5d05f5991 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -1217,14 +1217,11 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
         sharing the same x-coordinates can be (independently) fit with one
         call to `polyfit` by passing in for `y` a 2-D array that contains
         one data set per column.
-    deg : int or array_like
-        Degree of the fitting polynomial. If `deg` is a single integer
-        all terms up to and including the `deg`'th term are included.
-        `deg` may alternatively be a list or array specifying which
-        terms in the Legendre expansion to include in the fit.
-
-        .. versionchanged:: 1.11.0
-        `deg` may be a list specifying which terms to fit
+    deg : int or 1-D array_like
+        Degree(s) of the fitting polynomials. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included in the
+        fit. For Numpy versions >= 1.11 a list of integers specifying the
+        degrees of the terms to include may be used instead.
     rcond : float, optional
         Relative condition number of the fit.  Singular values smaller
         than `rcond`, relative to the largest singular value, will be
@@ -1340,11 +1337,11 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     """
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
-    deg = np.asarray([deg,], dtype=int).flatten()
+    deg = np.asarray(deg)
 
     # check arguments.
-    if deg.size < 1:
-        raise TypeError("expected deg to be one or more integers")
+    if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+        raise TypeError("deg must be an int or non-empty 1-D array of int")
     if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
@@ -1356,19 +1353,17 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
-    if deg.size == 1:
-        restricted_fit = False
-        lmax = deg[0]
+    if deg.ndim == 0:
+        lmax = deg
         order = lmax + 1
+        van = polyvander(x, lmax)
     else:
-        restricted_fit = True
-        lmax = deg.max()
-        order = deg.size
+        deg = np.sort(deg)
+        lmax = deg[-1]
+        order = len(deg)
+        van = polyvander(x, lmax)[:, deg]
 
     # set up the least squares matrices in transposed form
-    van = polyvander(x, lmax)
-    if restricted_fit:
-        van = van[:, deg]
     lhs = van.T
     rhs = y.T
     if w is not None:
@@ -1398,11 +1393,11 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     c = (c.T/scl).T
 
     # Expand c to include non-fitted coefficients which are set to zero
-    if restricted_fit:
+    if deg.ndim == 1:
         if c.ndim == 2:
-            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+            cc = np.zeros((lmax + 1, c.shape[1]), dtype=c.dtype)
         else:
-            cc = np.zeros(lmax+1, dtype=c.dtype)
+            cc = np.zeros(lmax + 1, dtype=c.dtype)
         cc[deg] = c
         c = cc