--- a/scripts/polynomial/polyfit.m	Wed Jun 24 07:42:37 2020 -0700
+++ b/scripts/polynomial/polyfit.m	Thu Jun 25 15:22:37 2020 -0700
@@ -27,25 +27,33 @@
 ## @deftypefn  {} {@var{p} =} polyfit (@var{x}, @var{y}, @var{n})
 ## @deftypefnx {} {[@var{p}, @var{s}] =} polyfit (@var{x}, @var{y}, @var{n})
 ## @deftypefnx {} {[@var{p}, @var{s}, @var{mu}] =} polyfit (@var{x}, @var{y}, @var{n})
-## Return the coefficients of a polynomial @var{p}(@var{x}) of degree
-## @var{n} that minimizes the least-squares-error of the fit to the points
-## @code{[@var{x}, @var{y}]}.
+## Return the coefficients of a polynomial @var{p}(@var{x}) of degree @var{n}
+## that minimizes the least-squares-error of the fit to the points
+## @code{[@var{x}(:), @var{y}(:)]}.
 ##
-## If @var{n} is a logical vector, it is used as a mask to selectively force
-## the corresponding polynomial coefficients to be used or ignored.
+## @var{n} is typically an integer @geq{} 0 specifying the degree of the
+## approximating polynomial.  If @var{n} is a logical vector, it is used as a
+## mask to selectively force the corresponding polynomial coefficients to be
+## used or ignored.
 ##
-## The polynomial coefficients are returned in a row vector.
+## The polynomial coefficients are returned in the row vector @var{p}.  The
+## output @var{p} may be directly used with @code{polyval} to estimate values
+## using the fitted polynomial.
 ##
 ## The optional output @var{s} is a structure containing the following fields:
 ##
 ## @table @samp
-## @item R
-## Triangular factor R from the QR@tie{}decomposition.
+##
+## @item yf
+## The values of the polynomial for each value of @var{x}.
 ##
 ## @item X
 ## The @nospell{Vandermonde} matrix used to compute the polynomial
 ## coefficients.
 ##
+## @item R
+## Triangular factor R from the QR@tie{}decomposition.
+##
 ## @item C
 ## The unscaled covariance matrix, formally equal to the inverse of
 ## @var{x'}*@var{x}, but computed in a way minimizing roundoff error
@@ -56,25 +64,37 @@
 ##
 ## @item normr
 ## The norm of the residuals.
-##
-## @item yf
-## The values of the polynomial for each value of @var{x}.
 ## @end table
 ##
-## The second output may be used by @code{polyval} to calculate the
-## statistical error limits of the predicted values.  In particular, the
-## standard deviation of @var{p} coefficients is given by
+## The second output may be used by @code{polyval} to calculate the statistical
+## error limits of the predicted values.  In particular, the standard deviation
+## of @var{p} coefficients is given by
 ##
-## @code{sqrt (diag (s.C)/s.df)*s.normr}.
+## @code{sqrt (diag (@var{s.C})/@var{s.df}) * @var{s.normr}}.
 ##
-## When the third output, @var{mu}, is present the coefficients, @var{p}, are
-## associated with a polynomial in
+## When the third output, @var{mu}, is present the original data is centered
+## and scaled which can improve the numerical stability of the fit.  The
+## coefficients @var{p} are associated with a polynomial in
 ##
 ## @code{@var{xhat} = (@var{x} - @var{mu}(1)) / @var{mu}(2)} @*
 ## where @var{mu}(1) = mean (@var{x}), and @var{mu}(2) = std (@var{x}).
 ##
-## This linear transformation of @var{x} improves the numerical stability of
-## the fit.
+## Example 1 : logical @var{n} and integer @var{n}
+##
+## @example
+## @group
+## f = @@(x) x.^2 + 5;   # data-generating function
+## x = 0:5;
+## y = f (x);
+## ## Fit data to polynomial A*x^3 + B*x^1
+## p = polyfit (x, y, logical ([1, 0, 1, 0]))
+## @result{} p = [ 0.0680, 0, 4.2444, 0 ]
+## ## Fit data to polynomial using all terms up to x^3
+## p = polyfit (x, y, 3)
+## @result{} p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]
+## @end group
+## @end example
+##
 ## @seealso{polyval, polyaffine, roots, vander, zscore}
 ## @end deftypefn
 
@@ -84,93 +104,118 @@
     print_usage ();
   endif
 
+  y_is_row_vector = isrow (y);
+
+  ## Reshape x & y into column vectors.
+  x = x(:);
+  y = y(:);
+
+  nx = numel (x);
+  ny = numel (y);
+  if (nx != ny)
+    error ("polyfit: X and Y must have the same number of points");
+  endif
+
   if (nargout > 2)
-    ## Normalized the x values.
+    ## Center and scale the x values.
     mu = [mean(x), std(x)];
     x = (x - mu(1)) / mu(2);
   endif
 
-  if (! size_equal (x, y))
-    error ("polyfit: X and Y must be vectors of the same size");
-  endif
-
+  ## n is the polynomial degree (an input, or deduced from the polymask size)
+  ## m is the effective number of coefficients.
   if (islogical (n))
-    polymask = n;
-    ## n is the polynomial degree as given the polymask size; m is the
-    ## effective number of used coefficients.
-    n = length (polymask) - 1; m = sum (polymask) - 1;
+    polymask = n(:).';          # force to row vector
+    n = numel (polymask) - 1;
+    m = sum (polymask) - 1;
+    pad_output = true;
   else
     if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == fix (n)))
       error ("polyfit: N must be a non-negative integer");
     endif
-    polymask = logical (ones (1, n+1)); m = n;
+    polymask = true (1, n+1);
+    m = n;
+    pad_output = false;
   endif
 
-  y_is_row_vector = (rows (y) == 1);
-
-  ## Reshape x & y into column vectors.
-  l = numel (x);
-  x = x(:);
-  y = y(:);
+  if (m >= nx)
+    warning ("polyfit: degree of polynomial N is >= number of data points; solution is not unique");
+    m = nx;
+    pad_output = true;
+    ## Keep the lowest m entries in polymask
+    idx = find (polymask);
+    idx((end-m+1):end) = [];  
+    polymask(idx) = false;
+  endif
 
   ## Construct the Vandermonde matrix.
-  v = vander (x, n+1);
+  X = vander (x, n+1);
+  v = X(:, polymask);
 
   ## Solve by QR decomposition.
-  [q, r, k] = qr (v(:, polymask), 0);
+  [q, r, k] = qr (v, 0);
   p = r \ (q' * y);
   p(k) = p;
 
-  if (n != m)
-    q = p; p = zeros (n+1, 1);
-    p(polymask) = q;
-  endif
-
-  if (nargout > 1)
+  if (isargout (2))
     yf = v*p;
-
     if (y_is_row_vector)
       s.yf = yf.';
     else
       s.yf = yf;
     endif
-    s.X = v;
+
+    s.X = X;
 
-    ## r.'*r is positive definite if X(:, polymask) is of full rank.
-    ## Invert it by cholinv to avoid taking the square root of squared
-    ## quantities.  If cholinv fails, then X(:, polymask) is rank deficient
-    ## and not invertible.
+    ## r.'*r is positive definite if matrix v is of full rank.  Invert it by
+    ## cholinv to avoid taking the square root of squared quantities.
+    ## If cholinv fails, then v is rank deficient and not invertible.
     try
       C = cholinv (r.'*r)(k, k);
     catch
-      C = NaN (m+1, m+1);
+      C = NaN (m, m);
     end_try_catch
 
-    if (n != m)
-      ## fill matrices if required
+    if (pad_output)
       s.X(:, ! polymask) = 0;
-      s.R = zeros (n+1, n+1); s.R(polymask, polymask) = r;
-      s.C = zeros (n+1, n+1); s.C(polymask, polymask) = C;
+      s.R = zeros (rows (r), n+1); s.R(:, polymask) = r;
+      s.C = zeros (rows (C), n+1); s.C(:, polymask) = C;
     else
       s.R = r;
       s.C = C;
     endif
-    s.df = l - m - 1;
+
+    s.df = max (0, nx - m - 1);
     s.normr = norm (yf - y);
   endif
 
-  ## Return a row vector.
-  p = p.';
+  if (pad_output)
+    ## Zero pad output
+    q = p;
+    p = zeros (n+1, 1);
+    p(polymask) = q;
+  endif
+  p = p.';  # Return a row vector.  
 
 endfunction
 
 
 %!shared x
 %! x = [-2, -1, 0, 1, 2];
-%!assert (polyfit (x, x.^2+x+1, 2), [1, 1, 1], sqrt (eps))
-%!assert (polyfit (x, x.^2+x+1, 3), [0, 1, 1, 1], sqrt (eps))
-%!fail ("polyfit (x, x.^2+x+1)")
-%!fail ("polyfit (x, x.^2+x+1, [])")
+
+%!assert (polyfit (x, 3*x.^2 + 2*x + 1, 2), [3, 2, 1], 10*eps)
+%!assert (polyfit (x, 3*x.^2 + 2*x + 1, logical ([1 1 1])), [3, 2, 1], 10*eps)
+%!assert (polyfit (x, x.^2 + 2*x + 3, 3), [0, 1, 2, 3], 10*eps)
+%!assert (polyfit (x, x.^2 + 2*x + 3, logical ([0 1 1 1])), [0 1 2 3], 10*eps)
+
+## Test logical input N
+%!test
+%! x = [0:5];
+%! y = 3*x.^3 + 2*x.^2 + 4;
+%! [p, s] = polyfit (x, y, logical ([1, 0, 1, 1]));
+%! assert (p(2), 0);
+%! assert (all (p([1, 3, 4])));
+%! assert (s.df, 3);
 
 ## Test difficult case where scaling is really needed.  This example
 ## demonstrates the rather poor result which occurs when the dependent
@@ -188,25 +233,29 @@
 %! assert (s2.normr < s1.normr);
 
 %!test
-%! x = 1:4;
-%! p0 = [1i, 0, 2i, 4];
-%! y0 = polyval (p0, x);
-%! p = polyfit (x, y0, numel (p0) - 1);
-%! assert (p, p0, 1000*eps);
-
-%!test
 %! warning ("off", "Octave:nearly-singular-matrix", "local");
 %! x = 1000 + (-5:5);
 %! xn = (x - mean (x)) / std (x);
 %! pn = ones (1,5);
 %! y = polyval (pn, xn);
-%! [p, s, mu] = polyfit (x, y, numel (pn) - 1);
-%! [p2, s2] = polyfit (x, y, numel (pn) - 1);
+%! n = numel (pn) - 1;
+%! [p, s, mu] = polyfit (x, y, n);
+%! [p2, s2] = polyfit (x, y, n);
 %! assert (p, pn, s.normr);
 %! assert (s.yf, y, s.normr);
 %! assert (mu, [mean(x), std(x)]);
 %! assert (s.normr/s2.normr < sqrt (eps));
 
+## Complex polynomials
+%!test
+%! x = 1:4;
+%! p0 = [1i, 0, 2i, 4];
+%! y = polyval (p0, x);
+%! n = numel (p0) - 1;
+%! p = polyfit (x, y, n);
+%! assert (p, p0, 1000*eps);
+
+## Matrix input
 %!test
 %! x = [1, 2, 3; 4, 5, 6];
 %! y = [0, 0, 1; 1, 0, 0];
@@ -214,4 +263,89 @@
 %! expected = [0, 1, -14, 65, -112, 60] / 12;
 %! assert (p, expected, sqrt (eps));
 
-%!error <vectors of the same size> polyfit ([1, 2; 3, 4], [1, 2, 3, 4], 2)
+## Orientation of output
+%!test
+%! x = 0:5;
+%! y = x.^4 + 2*x + 5;
+%! [p, s] = polyfit (x, y, 3);
+%! assert (isrow (s.yf));
+%! [p, s] = polyfit (x, y.', 3);
+%! assert (iscolumn (s.yf));
+
+## Insufficient data for fit
+%!test
+%! x = [1, 2];
+%! y = [3, 4];
+%! ## Disable warnings entirely because there is not a specific ID to disable.
+%! wstate = warning ();
+%! unwind_protect
+%!   warning ("off", "all");
+%!   p0 = polyfit (x, y, 4);
+%!   [p1, s, mu] = polyfit (x, y, 4);
+%! unwind_protect_cleanup
+%!   warning (wstate);
+%! end_unwind_protect
+%! assert (p0, [0, 0, 0, 1, 2], 10*eps);
+%! assert (p1, [0, 0, 0, sqrt(2)/2, 3.5], 10*eps);
+%! assert (size (s.X), [2, 5]);
+%! assert (s.X(:,1:3), zeros (2,3));
+%! assert (size (s.R), [2, 5]);
+%! assert (s.R(:,1:3), zeros (2,3));
+%! assert (size (s.C), [2, 5]);
+%! assert (s.C(:,1:3), zeros (2,3));
+%! assert (s.df, 0);
+%! assert (mu, [1.5, sqrt(2)/2]);
+
+%!test
+%! x = [1, 2, 3];
+%! y = 2*x + 1;
+%! ## Disable warnings entirely because there is not a specific ID to disable.
+%! wstate = warning ();
+%! unwind_protect
+%!   warning ("off", "all");
+%!   p0 = polyfit (x, y, logical ([1, 1, 1, 0 1]));
+%!   [p1, s, mu] = polyfit (x, y, logical ([1, 1, 1, 0 1]));
+%! unwind_protect_cleanup
+%!   warning (wstate);
+%! end_unwind_protect
+%! assert (p0, [0, -2/11, 12/11, 0, 23/11], 10*eps);
+%! assert (p1, [0, 2, 0, 0, 5], 10*eps);
+%! assert (size (s.X), [3, 5]);
+%! assert (s.X(:,[1,4]), zeros (3,2));
+%! assert (size (s.R), [3, 5]);
+%! assert (s.R(:,[1,4]), zeros (3,2));
+%! assert (size (s.C), [3, 5]);
+%! assert (s.C(:,[1,4]), zeros (3,2));
+%! assert (s.df, 0);
+%! assert (mu, [2, 1]);
+
+%!test <*57964>
+%! ## Disable warnings entirely because there is not a specific ID to disable.
+%! wstate = warning ();
+%! unwind_protect
+%!   warning ("off", "all");
+%!   [p, s] = polyfit ([1,2], [3,4], 2);
+%! unwind_protect_cleanup
+%!   warning (wstate);
+%! end_unwind_protect
+%! assert (size (p), [1, 3]);
+%! assert (size (s.X), [2, 3]);
+%! assert (s.X(:,1), [0; 0]);
+%! assert (size (s.R), [2, 3]);
+%! assert (s.R(:,1), [0; 0]);
+%! assert (size (s.C), [2, 3]);
+%! assert (s.C(:,1), [0; 0]);
+
+## Test input validation
+%!error polyfit ()
+%!error polyfit (1)
+%!error polyfit (1,2)
+%!error polyfit (1,2,3,4,5)
+%!error <X and Y must have the same number of points> polyfit ([1, 2], 1, 1)
+%!error <X and Y must have the same number of points> polyfit (1, [1, 2], 1)
+%!error <N must be a non-negative integer> polyfit (1, 2, [1,2])
+%!error <N must be a non-negative integer> polyfit (1, 2, -1)
+%!error <N must be a non-negative integer> polyfit (1, 2, Inf)
+%!error <N must be a non-negative integer> polyfit (1, 2, 1.5)
+%!test <*57964>
+%! fail ("p = polyfit ([1,2], [3,4], 4)", "warning", "solution is not unique");