# HG changeset patch
# User Pascal Dupuis <Pascal.Dupuis@uclouvain.be>
# Date 1328450865 -3600
# Node ID c87f927d047da0e96535c14ae00d76bf97009d90
# Parent  7dd6ac033e69c20e6f993d0f5cd34da720e84b4c
polyfit.m: add the ability to specify the polynomial template.

diff --git a/scripts/polynomial/polyfit.m b/scripts/polynomial/polyfit.m
--- a/scripts/polynomial/polyfit.m
+++ b/scripts/polynomial/polyfit.m
@@ -22,7 +22,9 @@
 ## @deftypefnx {Function File} {[@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}]}.
+## @code{[@var{x}, @var{y}]}. If @var{n} is a logical vector, it is used
+## as a mask to selectivelly force the corresponding polynomial
+## coefficients to be used or ignored.
 ##
 ## The polynomial coefficients are returned in a row vector.
 ##
@@ -35,6 +37,11 @@
 ## @item X
 ## The Vandermonde matrix used to compute the polynomial coefficients.
 ##
+## @item C
+## The unscaled covariance matrix, formally equal to the inverse of
+## @var{x'}*@var{x}, but computed in a way minimizing roundoff error
+## propagation.
+## 
 ## @item df
 ## The degrees of freedom.
 ##
@@ -46,7 +53,9 @@
 ## @end table
 ##
 ## The second output may be used by @code{polyval} to calculate the
-## statistical error limits of the predicted values.
+## 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}.
 ##
 ## When the third output, @var{mu}, is present the
 ## coefficients, @var{p}, are associated with a polynomial in
@@ -60,6 +69,8 @@
 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
 ## Created: 13 December 1994
 ## Adapted-By: jwe
+## Modified on 20120204 by P. Dupuis; added the ability to specify a
+## polynomial mask instead of a polynomial degree.
 
 function [p, s, mu] = polyfit (x, y, n)
 
@@ -77,8 +88,16 @@
     error ("polyfit: X and Y must be vectors of the same size");
   endif
 
-  if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == fix (n)))
-    error ("polyfit: N must be a non-negative integer");
+  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;
+  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;
   endif
 
   y_is_row_vector = (rows (y) == 1);
@@ -92,10 +111,15 @@
   v = vander (x, n+1);
 
   ## Solve by QR decomposition.
-  [q, r, k] = qr (v, 0);
+  [q, r, k] = qr (v(:, polymask), 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)
     yf = v*p;
 
@@ -104,10 +128,28 @@
     else
       s.yf = yf;
     endif
+    s.X = v; 
 
-    s.R = r;
-    s.X = v;
-    s.df = l - n - 1;
+    ## 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.
+    try
+      C = cholinv (r.'*r)(k, k);
+    catch
+      C = NaN * ones (m+1, m+1);
+    end_try_catch
+
+    if (n ~= m)
+      ## fill matrices if required
+      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;
+    else
+      s.R = r; 
+      s.C = C;
+    endif
+    s.df = l - m - 1;
     s.normr = norm (yf - y);
   endif