# HG changeset patch
# User jwe
# Date 1148021611 0
# Node ID 448f9982e7fbed1f40b1ee1afecf574159c0a62c
# Parent  080c08b192d8d0c62bcc8c069843185b5372b1de
[project @ 2006-05-19 06:53:31 by jwe]

diff -r 080c08b192d8 -r 448f9982e7fb scripts/ChangeLog
--- a/scripts/ChangeLog	Fri May 19 05:32:19 2006 +0000
+++ b/scripts/ChangeLog	Fri May 19 06:53:31 2006 +0000
@@ -1,3 +1,8 @@
+2006-05-19  David Bateman  <dbateman@free.fr>
+
+	* polynomial/unmkpp.m, polynomial/mkpp.m, polynomial/spline.m,
+	polynomial/ppval.m: New files from Octave Forge.
+
 2006-05-17  John W. Eaton  <jwe@octave.org>
 
 	* set/intersection.m: Delete
diff -r 080c08b192d8 -r 448f9982e7fb scripts/polynomial/mkpp.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/mkpp.m	Fri May 19 06:53:31 2006 +0000
@@ -0,0 +1,67 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pp} = } mkpp (@var{x}, @var{p})
+## @deftypefnx {Function File} {@var{pp} = } mkpp (@var{x}, @var{p}, @var{d})
+## 
+## Construct a piece-wise polynomial structure from sample points
+## @var{x} and coefficients @var{p}.  The ith row of @var{p},
+## @code{@var{p} (@var{i},:)}, contains the coefficients for the polynomial
+## over the @var{i}-th interval, ordered from highest to 
+## lowest. There must be one row for each interval in @var{x}, so 
+## @code{rows (@var{p}) == length (@var{x}) - 1}.  
+##
+## You can concatenate multiple polynomials of the same order over the 
+## same set of intervals using @code{@var{p} = [ @var{p1}; @var{p2}; 
+## @dots{}; @var{pd} ]}.  In this case, @code{rows (@var{p}) == @var{d} 
+## * (length (@var{x}) - 1)}.
+##
+## @var{d} specifies the shape of the matrix @var{p} for all except the
+## last dimension. If @var{d} is not specified it will be computed as
+## @code{round (rows (@var{p}) / (length (@var{x}) - 1)) instead.
+##
+## @seealso{unmkpp, ppval, spline}
+## @end deftypefn
+
+function pp = mkpp (x, P, d)
+  if (nargin < 2 || nargin > 3)
+    usage ("pp = mkpp(x,P,d)");
+  endif
+  pp.x = x(:);
+  pp.P = P;
+  pp.n = length (x) - 1;
+  pp.k = columns (P);
+  if (nargin < 3)
+    d = round (rows (P) / pp.n); 
+  endif
+  pp.d = d;
+  if (pp.n*d != rows (P))
+    error ("mkpp: num intervals in x doesn't match num polynomials in P");
+  endif
+endfunction
+
+%!demo # linear interpolation
+%! x=linspace(0,pi,5)'; 
+%! t=[sin(x),cos(x)];
+%! m=diff(t)./(x(2)-x(1)); 
+%! b=t(1:4,:);
+%! pp = mkpp(x, [m(:),b(:)]);
+%! xi=linspace(0,pi,50);
+%! plot(x,t,"x;control;",xi,ppval(pp,xi),";interp;");
diff -r 080c08b192d8 -r 448f9982e7fb scripts/polynomial/ppval.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/ppval.m	Fri May 19 06:53:31 2006 +0000
@@ -0,0 +1,58 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp, @var{xi})
+##
+## Evaluate piece-wise polynomial @var{pp} at the points @var{xi}.  
+## If @code{@var{pp}.d} is a scalar greater than 1, or an array,
+## then the returned value @var{yi} will be an array that is 
+## @code{d1, d1, @dots{}, dk, length (@var{xi})]}.
+##
+## @seealso{mkpp, unmkpp, spline}
+## @end deftypefn 
+
+function yi = ppval (pp, xi)
+
+  if (nargin != 2)
+    usage ("yi = ppval(pp, xi)")
+  endif
+  if (! isstruct (pp))
+    error ("ppval: expects a pp structure");
+  endif
+  if (isempty (xi))
+    yi = [];
+  else
+    transposed = (columns (xi) == 1);
+    xi = xi(:);
+    xn = length (xi);
+    idx = lookup (pp.x(2:pp.n), xi) + 1;
+    dx = (xi - pp.x(idx))';
+    dx = reshape (dx(ones(1,prod(pp.d)),:),[pp.d,xn]);
+    c = reshape (pp.P(:,1), pp.n, prod (pp.d));
+    yi = reshape (c(idx,:)', [pp.d, xn]);
+    for i  = 2 : pp.k;
+      c = reshape (pp.P(:,i), pp.n, prod (pp.d));
+      yi = yi .* dx + reshape (c(idx,:)', [pp.d, xn]);
+    endfor
+    if (transposed && isscalar (pp.d) && pp.d == 1)
+      yi = yi.';
+    endif
+  endif
+endfunction
diff -r 080c08b192d8 -r 448f9982e7fb scripts/polynomial/spline.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/spline.m	Fri May 19 06:53:31 2006 +0000
@@ -0,0 +1,224 @@
+## Copyright (C) 2000,2001  Kai Habel
+## Copyright (C) 2006 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pp} = } spline (@var{x}, @var{y})
+## @deftypefnx {Function File} {@var{yi} = } spline (@var{x}, @var{y}, @var{xi})
+##
+## Returns the cubic spline interpolation of @var{y} at the point
+## @var{x}. Called with two arguments the piece-wse polynomial @var{pp}
+## that may later be used with @code{ppval} to evaluate the polynomial
+## at specific points.
+##
+## The variable @var{x} must be a vector of length @var{n}, and @var{y}
+## can be either a vector of array. In the case where @var{y} is a
+## vector, it can have a length of either @var{n} or @code{@var{n} + 2}.
+## If the length of @var{y} is @var{n}, then the 'not-a-knot' end
+## condition is used. If the length of @var{y} is @code{@var{n} + 2},
+## then the first and last values of the vector @var{y} are the first
+## derivative of the cubic spline at the end-points.
+##
+## If @var{y} is an array, then the size of @var{y} must have the form
+## @iftex
+## @tex
+## $$[s_1, s_2, \cdots, s_k, n]$$
+## @end tex
+## @end iftex
+## @ifinfo
+## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
+## @end ifinfo
+## or
+## @iftex
+## @tex
+## $$[s_1, s_2, \cdots, s_k, n + 2]$$.
+## @end tex
+## @end iftex
+## @ifinfo
+## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n} + 2]}.
+## @end ifinfo
+## The array is then reshaped internally to a matrix where to leading
+## dimension is given by 
+## @iftex
+## @tex
+## $$s_1 s_2 \cdots s_k$$
+## @end tex
+## @end iftex
+## @ifinfo
+## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
+## @end ifinfo
+## and each row this matrix is then treated seperately. Note that this
+## is exactly the opposite treatment than @code{interp1} and is done
+## for compatiability.
+##
+## Called with a third input argument, @code{spline} evaluates the 
+## piece-wise spline at the points @var{xi}. There is an equivalence
+## between @code{ppval (spline (@var{x}, @var{y}), @var{xi})} and
+## @code{spline (@var{x}, @var{y}, @var{xi})}.
+##
+## @seealso{ppval, mkpp, unmkpp}
+## @end deftypefn
+
+## This code is based on csape.m from octave-forge, but has been
+## modified to use the sparse solver code in octave that itself allows
+## special casing of tri-diagonal matrices, modified for NDArrays and
+## for the treatment of vectors y 2 elements longer than x as complete
+## splines.
+
+function ret = spline (x, y, xi)
+
+  x = x(:);
+  n = length (x);
+  if (n < 3)
+    error ("spline: requires at least 3 points"); 
+  endif
+
+  ## Check the size and shape of y
+  ndy = ndims (y);
+  szy = size (y);
+  if (ndy == 2 && (szy(1) == 1 || szy(2) == 1))
+    if (szy(1) == 1)
+      a = y';
+    else
+      a = y;
+      szy = fliplr (szy);
+    endif
+  else
+    a = reshape (y, [prod(szy(1:end-1)), szy(end)])';
+  endif
+  complete = false;
+  if (size (a, 1) == n + 2)
+    complete = true;
+    dfs = a(1,:);
+    dfe = a(end,:);
+    a = a(2:end-1,:);
+  endif
+
+  b = c = zeros (size (a));
+  h = diff (x);
+  idx = ones (columns (a),1);
+
+  if (complete)
+
+    if (n == 3)
+      dg = 1.5 * h(1) - 0.5 * h(2);
+      c(2:n - 1,:) = 1/dg(1);
+    else
+      dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
+      dg(1) = dg(1) - 0.5 * h(1);
+      dg(n - 2) = dg(n-2) - 0.5 * h(n - 1);
+
+      e = h(2:n - 2);
+
+      size(a)
+      size(h)
+      n
+
+      g = 3 * diff (a(2:n,:)) ./ h(2:n - 1,idx)\
+        - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2,idx);
+      g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) \
+          - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs);
+      g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - dfe)\
+          - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2);
+
+      c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g;
+    endif
+
+    c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
+	      - c(2,:) * h(1)) / (2 * h(1)); 
+    c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * dfe
+		+ c(n - 1,:) * h(n - 1)) / (2 * h(n - 1));
+    b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)\
+      - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
+    d = diff (c) ./ (3 * h(1:n - 1, idx));
+
+  else
+
+    g = zeros(n - 2,columns(a));
+    g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)\
+				  - h(2) / h(1) * (a(2,:) - a(1,:)));
+    g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *\
+	(h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\
+	 (a(n - 1,:) - a(n - 2,:)));
+
+    if (n > 4)
+
+      g(2:n - 3,:) = 3 * diff (a(3:n - 1,:)) ./ h(3:n - 2,idx)\
+          - 3 * diff (a(2:n - 2,:)) ./ h(2:n - 3,idx);
+
+      dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
+      dg(1) = dg(1) - h(1);
+      dg(n - 2) = dg(n-2) - h(n - 1);
+
+      ldg = udg = h(2:n - 2);
+      udg(1) = udg(1) - h(1);
+      ldg(n - 3) = ldg(n-3) - h(n - 1);
+      c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
+
+    elseif (n == 4)
+
+      dg = [h(1) + 2 * h(2), 2 * h(2) + h(3)];
+      ldg = h(2) - h(3);
+      udg = h(2) - h(1);
+      c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
+      
+    else # n == 3
+	    
+      dg= [h(1) + 2 * h(2)];
+      c(2:n - 1,:) = g/dg(1);
+
+    endif
+
+    c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
+    c(n,:) = c(n - 1,:) + h(n - 1) / h(n - 2) * (c(n - 1,:) - c(n - 2,:));
+    b = diff (a) ./ h(1:n - 1, idx)\
+	- h(1:n - 1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
+    d = diff (c) ./ (3 * h(1:n - 1, idx));
+
+  endif
+
+  d = d(1:n-1,:); c=c(1:n-1,:); b=b(1:n-1,:); a=a(1:n-1,:);
+  coeffs = [d(:), c(:), b(:), a(:)];
+  ret = mkpp (x, coeffs, szy(1:end-1));
+
+  if (nargin == 3)
+    ret = ppval (ret, xi);
+  endif
+
+endfunction
+
+%!demo
+%! x = 0:10; y = sin(x);
+%! xspline = 0:0.1:10; yspline = spline(x,y,xspline);
+%! title("spline fit to points from sin(x)");
+%! plot(xspline,sin(xspline),";original;",...
+%!      xspline,yspline,"-;interpolation;",...
+%!      x,y,"+;interpolation points;");
+%! %--------------------------------------------------------
+%! % confirm that interpolated function matches the original
+
+%!shared x,y
+%! x = [0:10]; y = sin(x);
+%!assert (spline(x,y,x), y);
+%!assert (spline(x,y,x'), y');
+%!assert (spline(x',y',x'), y');
+%!assert (spline(x',y',x), y);
+%!assert (isempty(spline(x',y',[])));
+%!assert (isempty(spline(x,y,[])));
+%!assert (spline(x,[y;y],x), [spline(x,y,x);spline(x,y,x)])
diff -r 080c08b192d8 -r 448f9982e7fb scripts/polynomial/unmkpp.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/unmkpp.m	Fri May 19 06:53:31 2006 +0000
@@ -0,0 +1,61 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{x}, @var{p}, @var{n}, @var{k}, @var{d}] =} unmkpp (@var{pp})
+##
+## Extract the components of a piece-wise polynomial structure @var{pp}.
+## These are as follows:
+##
+## @table @asis
+## @item @var{x}
+## Samples points.
+## @item @var{p}
+## Polynomial coefficients for points in sample interval. @code{@var{p}
+## (@var{i}, :)} contains the coefficients for the polynomial over
+## interval @var{i} ordered from highest to lowest. If @code{@var{d} >
+## 1}, @code{@var{p} (@var{r}, @var{i}, :)} contains the coeffients for 
+## the r-th polynomial defined on interval @var{i}. However, this is 
+## stored as a 2-D array such that @code{@var{c} = reshape (@var{p} (:,
+## @var{j}), @var{n}, @var{d})} gives @code{@var{c} (@var{i},  @var{r})}
+## is the j-th coefficient of the r-th polynomial over the i-th interval.
+## @item @var{n}
+## Number of polynomial pieces.
+## @item @var{k}
+## Order of the polynomial plus 1.
+## @item @var{d}
+## Number of polynomials defined for each interval.
+## @end table
+##
+## @seealso{mkpp, ppval, spline}
+## @end deftypefn
+
+function [x, P, n, k, d] = unmkpp (pp)
+  if (nargin == 0)
+    usage ("[x, P, n, k, d] = unmkpp(pp)")
+  endif
+  if (! isstruct (pp))
+    error ("unmkpp: expecting piecewise polynomial structure");
+  endif
+  x = pp.x;
+  P = pp.P;
+  n = pp.n;
+  k = pp.k;
+  d = pp.d;
+endfunction