Mercurial > octave-nkf
changeset 14509:a88f8e4fae56
New Function, splinefit.m
* __splinefit__.m: New private file. Jonas Lundgren's splinefit.m with BSD
License. Jonas emailed this version to Octave's developers.
* splinefit.m: New File. Piece-wise polynomial fit. This is a wrapper for
__splinefit__.m. The wrapper allows for Octave's tex-info documentation,
demos, and tests to be added. In addition the input syntax has been sligtly
modified to allow new options to be added without breaking compatiblity.
* doc/splineimages.m: New file to produce splineimages<#> for the docs.
* doc/images: Include splineimages.m and the figues for the docs.
* scripts/polynomial/module.mk: Add new files.
* doc/interpreter/poly.txi: Minimal description of splinefit.
| author | Ben Abbott <bpabbott@mac.com> |
|---|---|
| date | Thu, 29 Mar 2012 19:13:21 -0400 |
| parents | 0901f926ed50 |
| children | c20de00a66a9 |
| files | doc/interpreter/images doc/interpreter/poly.txi doc/interpreter/splineimages.m scripts/polynomial/module.mk scripts/polynomial/private/__splinefit__.m scripts/polynomial/splinefit.m |
| diffstat | 6 files changed, 1272 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/interpreter/images Wed Mar 28 23:09:43 2012 -0400 +++ b/doc/interpreter/images Thu Mar 29 19:13:21 2012 -0400 @@ -2,3 +2,4 @@ interpimages.m interpft interpn interpderiv1 interpderiv2 plotimages.m plot hist errorbar polar mesh plot3 extended sparseimages.m gplot grid spmatrix spchol spcholperm +splineimages.m splinefit1 splinefit2 splinefit3 splinefit4 splinefit6
--- a/doc/interpreter/poly.txi Wed Mar 28 23:09:43 2012 -0400 +++ b/doc/interpreter/poly.txi Thu Mar 29 19:13:21 2012 -0400 @@ -132,18 +132,228 @@ Octave comes with good support for various kinds of interpolation, most of which are described in @ref{Interpolation}. One simple alternative to the functions described in the aforementioned chapter, is to fit -a single polynomial to some given data points. To avoid a highly -fluctuating polynomial, one most often wants to fit a low-order polynomial -to data. This usually means that it is necessary to fit the polynomial -in a least-squares sense, which just is what the @code{polyfit} function does. +a single polynomial, or a piecewise polynomial (spline) to some given +data points. To avoid a highly fluctuating polynomial, one most often +wants to fit a low-order polynomial to data. This usually means that it +is necessary to fit the polynomial in a least-squares sense, which just +is what the @code{polyfit} function does. @DOCSTRING(polyfit) In situations where a single polynomial isn't good enough, a solution -is to use several polynomials pieced together. The function @code{mkpp} -creates a piecewise polynomial, @code{ppval} evaluates the function -created by @code{mkpp}, and @code{unmkpp} returns detailed information -about the function. +is to use several polynomials pieced together. The function +@code{splinefit} fits a peicewise polynomial (spline) to a set of +data. + +@DOCSTRING(splinefit) + +The number of @var{breaks} (or knots) used to construct the piecewise +polynomial is a significant factor in suppressing the noise present in +the input data, @var{x} and @var{y}. This is demostrated by the example +below. + +@example +@group +x = 2 * pi * rand (1, 200); +y = sin (x) + sin (2 * x) + 0.2 * randn (size (x)); +## Uniform breaks +breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces +pp1 = splinefit (x, y, breaks); +## Breaks interpolated from data +pp2 = splinefit (x, y, 10); % 11 breaks, 10 pieces +## Plot +xx = linspace (0, 2 * pi, 400); +y1 = ppval (pp1, xx); +y2 = ppval (pp2, xx); +plot (x, y, ".", xx, [y1; y2]) +axis tight +ylim auto +legend (@{"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"@}) +@end group +@end example + +@ifnotinfo +@noindent +The result of which can be seen in @ref{fig:splinefit1}. + +@float Figure,fig:splinefit1 +@center @image{splinefit1,4in} +@caption{Comparison of a fitting a piecewise polynomial with 41 breaks to one +with 11 breaks. The fit with the large number of breaks exhibits a fast ripple +that is not present in the underlying function.} +@end float +@end ifnotinfo + +The piece-wise polynomial fit provided by @code{splinefit} provides +continuous derivatives up to the @var{order} of the fit. This can +be demonstrated by the code + +@example +@group +## Data (200 points) +x = 2 * pi * rand (1, 200); +y = sin (x) + sin (2 * x) + 0.1 * randn (size (x)); +## Piecewise constant +pp1 = splinefit (x, y, 8, "order", 0); +## Piecewise linear +pp2 = splinefit (x, y, 8, "order", 1); +## Piecewise quadratic +pp3 = splinefit (x, y, 8, "order", 2); +## Piecewise cubic +pp4 = splinefit (x, y, 8, "order", 3); +## Piecewise quartic +pp5 = splinefit (x, y, 8, "order", 4); +## Plot +xx = linspace (0, 2 * pi, 400); +y1 = ppval (pp1, xx); +y2 = ppval (pp2, xx); +y3 = ppval (pp3, xx); +y4 = ppval (pp4, xx); +y5 = ppval (pp5, xx); +plot (x, y, ".", xx, [y1; y2; y3; y4; y5]) +axis tight +ylim auto +legend (@{"data", "order 0", "order 1", "order 2", "order 3", "order 4"@}) +@end group +@end example + +@ifnotinfo +@noindent +The result of which can be seen in @ref{fig:splinefit2}. + +@float Figure,fig:splinefit2 +@center @image{splinefit2,4in} +@caption{Comparison of a piecewise constant, linear, quadratic, cubic, and +quartic polynomials with 8 breaks to noisey data. The higher order solutions +more accurately represent the underlying function, but come with the +expense of computational complexity.} +@end float +@end ifnotinfo + +When the underlying function to provide a fit to is periodice, @code{splitfit} +is able to apply the boundary conditions needed to manifest a periodic fit. +This is demonstrated by the code below. + +@example +@group +## Data (100 points) +x = 2 * pi * [0, (rand (1, 98)), 1]; +y = sin (x) - cos (2 * x) + 0.2 * randn (size (x)); +## No constraints +pp1 = splinefit (x, y, 10, "order", 5); +## Periodic boundaries +pp2 = splinefit (x, y, 10, "order", 5, "periodic", true); +## Plot +xx = linspace (0, 2 * pi, 400); +y1 = ppval (pp1, xx); +y2 = ppval (pp2, xx); +plot (x, y, ".", xx, [y1; y2]) +axis tight +ylim auto +legend (@{"data", "no constraints", "periodic"@}) +@end group +@end example + +@ifnotinfo +@noindent +The result of which can be seen in @ref{fig:splinefit3}. + +@float Figure,fig:splinefit3 +@center @image{splinefit3,4in} +@caption{Comparison of piecewise cubic fits to a noisy periodic function with, +and without, periodic boundary conditions.} +@end float +@end ifnotinfo + +More complex constaints may be added as well. For example, the code below +illustrates a periodic fit with values that have been clamped end point values +and a second periodic fit wihh hinged end point values. + +@example +@group +## Data (200 points) +x = 2 * pi * rand (1, 200); +y = sin (2 * x) + 0.1 * randn (size (x)); +## Breaks +breaks = linspace (0, 2 * pi, 10); +## Clamped endpoints, y = y" = 0 +xc = [0, 0, 2*pi, 2*pi]; +cc = [(eye (2)), (eye (2))]; +con = struct ("xc", xc, "cc", cc); +pp1 = splinefit (x, y, breaks, "constraints", con); +## Hinged periodic endpoints, y = 0 +con = struct ("xc", 0); +pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true); +## Plot +xx = linspace (0, 2 * pi, 400); +y1 = ppval (pp1, xx); +y2 = ppval (pp2, xx); +plot (x, y, ".", xx, [y1; y2]) +axis tight +ylim auto +legend (@{"data", "clamped", "hinged periodic"@}) +@end group +@end example + +@ifnotinfo +@noindent +The result of which can be seen in @ref{fig:splinefit4}. + +@float Figure,fig:splinefit4 +@center @image{splinefit4,4in} +@caption{Comparison of two periodic piecewise cubic fits to a noisy periodic +signal. One fit has its end points clamped and the second has its end points +hinged.} +@end float +@end ifnotinfo + +The @code{splinefit} function also provides the convenience of a @var{robust} +fitting, where the effect of outlying data is reduced. In the example below, +three different fits are provided. Two with differing levels of outlier +suppressin and a third illustrating the non-robust solution. + +@example +@group +## Data +x = linspace (0, 2*pi, 200); +y = sin (x) + sin (2 * x) + 0.05 * randn (size (x)); +## Add outliers +x = [x, linspace(0,2*pi,60)]; +y = [y, -ones(1,60)]; +## Fit splines with hinged conditions +con = struct ("xc", [0, 2*pi]); +## Robust fitting, beta = 0.25 +pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25); +## Robust fitting, beta = 0.75 +pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75); +## No robust fitting +pp3 = splinefit (x, y, 8, "constraints", con); +## Plot +xx = linspace (0, 2*pi, 400); +y1 = ppval (pp1, xx); +y2 = ppval (pp2, xx); +y3 = ppval (pp3, xx); +plot (x, y, ".", xx, [y1; y2; y3]) +legend (@{"data with outliers","robust, beta = 0.25", ... + "robust, beta = 0.75", "no robust fitting"@}) +axis tight +ylim auto +@end group +@end example + +@ifnotinfo +@noindent +The result of which can be seen in @ref{fig:splinefit6}. + +@float Figure,fig:splinefit6 +@center @image{splinefit6,4in} +@caption{Comparison of two differnet levels of robust fitting (@var{beta} = 0.25 and 0.75) to noisy data combined with outlying data. A standard fit is also included.} +@end float +@end ifnotinfo + +The function, @code{ppval}, evaluates the piecewise polyomials, created +by @code{mkpp} or other means, and @code{unmkpp} returns detailed +information about the piecewise polynomial. The following example shows how to combine two linear functions and a quadratic into one function. Each of these functions is expressed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/doc/interpreter/splineimages.m Thu Mar 29 19:13:21 2012 -0400 @@ -0,0 +1,192 @@ +## Copyright (C) 2012 Ben Abbott, Jonas Lundgren +## +## 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 3 of the License, 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, see +## <http://www.gnu.org/licenses/>. + +function splineimages (nm, typ) + graphics_toolkit ("gnuplot"); + set_print_size (); + hide_output (); + if (strcmp (typ, "png")) + set (0, "defaulttextfontname", "*"); + endif + if (strcmp (typ, "eps")) + d_typ = "-depsc2"; + else + d_typ = cstrcat ("-d", typ); + endif + + if (strcmp (typ, "txt")) + image_as_txt (nm); + elseif (strcmp (nm, "splinefit1")) ## Breaks and Pieces + x = 2 * pi * rand (1, 200); + y = sin (x) + sin (2 * x) + 0.2 * randn (size (x)); + ## Uniform breaks + breaks = linspace (0, 2 * pi, 41); ## 41 breaks, 40 pieces + pp1 = splinefit (x, y, breaks); + ## Breaks interpolated from data + pp2 = splinefit (x, y, 10); ## 11 breaks, 10 pieces + ## Plot + xx = linspace (0, 2 * pi, 400); + y1 = ppval (pp1, xx); + y2 = ppval (pp2, xx); + plot (x, y, ".", xx, [y1; y2]) + axis tight + ylim ([-2.5 2.5]) + legend ("data", "41 breaks, 40 pieces", "11 breaks, 10 pieces") + print (cstrcat (nm, ".", typ), d_typ) + elseif (strcmp (nm, "splinefit2")) ## Spline orders + ## Data (200 points) + x = 2 * pi * rand (1, 200); + y = sin (x) + sin (2 * x) + 0.1 * randn (size (x)); + ## Splines + pp1 = splinefit (x, y, 8, "order", 0); ## Piecewise constant + pp2 = splinefit (x, y, 8, "order", 1); ## Piecewise linear + pp3 = splinefit (x, y, 8, "order", 2); ## Piecewise quadratic + pp4 = splinefit (x, y, 8, "order", 3); ## Piecewise cubic + pp5 = splinefit (x, y, 8, "order", 4); ## Etc. + ## Plot + xx = linspace (0, 2 * pi, 400); + y1 = ppval (pp1, xx); + y2 = ppval (pp2, xx); + y3 = ppval (pp3, xx); + y4 = ppval (pp4, xx); + y5 = ppval (pp5, xx); + plot (x, y, ".", xx, [y1; y2; y3; y4; y5]) + axis tight + ylim ([-2.5 2.5]) + legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"}) + print (cstrcat (nm, ".", typ), d_typ) + elseif (strcmp (nm, "splinefit3")) + ## Data (100 points) + x = 2 * pi * [0, (rand (1, 98)), 1]; + y = sin (x) - cos (2 * x) + 0.2 * randn (size (x)); + ## No constraints + pp1 = splinefit (x, y, 10, "order", 5); + ## Periodic boundaries + pp2 = splinefit (x, y, 10, "order", 5, "periodic", true); + ## Plot + xx = linspace (0, 2 * pi, 400); + y1 = ppval (pp1, xx); + y2 = ppval (pp2, xx); + plot (x, y, ".", xx, [y1; y2]) + axis tight + ylim ([-2 3]) + legend ({"data", "no constraints", "periodic"}) + print (cstrcat (nm, ".", typ), d_typ) + elseif (strcmp (nm, "splinefit4")) + ## Data (200 points) + x = 2 * pi * rand (1, 200); + y = sin (2 * x) + 0.1 * randn (size (x)); + ## Breaks + breaks = linspace (0, 2 * pi, 10); + ## Clamped endpoints, y = y" = 0 + xc = [0, 0, 2*pi, 2*pi]; + cc = [(eye (2)), (eye (2))]; + con = struct ("xc", xc, "cc", cc); + pp1 = splinefit (x, y, breaks, "constraints", con); + ## Hinged periodic endpoints, y = 0 + con = struct ("xc", 0); + pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true); + ## Plot + xx = linspace (0, 2 * pi, 400); + y1 = ppval (pp1, xx); + y2 = ppval (pp2, xx); + plot (x, y, ".", xx, [y1; y2]) + axis tight + ylim ([-1.5 1.5]) + legend({"data", "clamped", "hinged periodic"}) + print (cstrcat (nm, ".", typ), d_typ) + elseif (strcmp (nm, "splinefit5")) + ## Truncated data + x = [0, 1, 2, 4, 8, 16, 24, 40, 56, 72, 80] / 80; + y = [0, 28, 39, 53, 70, 86, 90, 79, 55, 22, 2] / 1000; + xy = [x; y]; + ## Curve length parameter + ds = sqrt (diff (x).^2 + diff (y).^2); + s = [0, cumsum(ds)]; + ## Constraints at s = 0: (x,y) = (0,0), (dx/ds,dy/ds) = (0,1) + con = struct ("xc", [0 0], "yc", [0 0; 0 1], "cc", eye (2)); + ## Fit a spline with 4 pieces + pp = splinefit (s, xy, 4, "constraints", con); + ## Plot + ss = linspace (0, s(end), 400); + xyfit = ppval (pp, ss); + xyb = ppval(pp, pp.breaks); + plot (x, y, ".", xyfit(1,:), xyfit(2,:), "r", xyb(1,:), xyb(2,:), "ro") + legend ({"data", "spline", "breaks"}) + axis tight + ylim ([0 0.1]) + print (cstrcat (nm, ".", typ), d_typ) + elseif (strcmp (nm, "splinefit6")) + ## Data + x = linspace (0, 2*pi, 200); + y = sin (x) + sin (2 * x) + 0.05 * randn (size (x)); + ## Add outliers + x = [x, linspace(0,2*pi,60)]; + y = [y, -ones(1,60)]; + ## Fit splines with hinged conditions + con = struct ("xc", [0, 2*pi]); + pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25); ## Robust fitting + pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75); ## Robust fitting + pp3 = splinefit (x, y, 8, "constraints", con); ## No robust fitting + ## Plot + xx = linspace (0, 2*pi, 400); + y1 = ppval (pp1, xx); + y2 = ppval (pp2, xx); + y3 = ppval (pp3, xx); + plot (x, y, ".", xx, [y1; y2; y3]) + legend({"data with outliers","robust, beta = 0.25", ... + "robust, beta = 0.75", "no robust fitting"}) + axis tight + ylim ([-2 2]) + print (cstrcat (nm, ".", typ), d_typ) + endif + hide_output (); +endfunction + +function set_print_size () + image_size = [5.0, 3.5]; # in inches, 16:9 format + border = 0; # For postscript use 50/72 + set (0, "defaultfigurepapertype", "<custom>"); + set (0, "defaultfigurepaperorientation", "landscape"); + set (0, "defaultfigurepapersize", image_size + 2*border); + set (0, "defaultfigurepaperposition", [border, border, image_size]); +endfunction + +## Use this function before plotting commands and after every call to +## print since print() resets output to stdout (unfortunately, gnpulot +## can't pop output as it can the terminal type). +function hide_output () + f = figure (1); + set (f, "visible", "off"); +endfunction + +## generate something for the texinfo @image command to process +function image_as_txt(nm) + fid = fopen (sprintf ("##s.txt", nm), "wt"); + fputs (fid, "\n"); + fputs (fid, "+---------------------------------+\n"); + fputs (fid, "| Image unavailable in text mode. |\n"); + fputs (fid, "+---------------------------------+\n"); + fclose (fid); +endfunction + +%!demo +%! for s = 1:6 +%! splineimages (sprintf ("splinefit##d", s), "pdf") +%! endfor +
--- a/scripts/polynomial/module.mk Wed Mar 28 23:09:43 2012 -0400 +++ b/scripts/polynomial/module.mk Thu Mar 29 19:13:21 2012 -0400 @@ -1,5 +1,8 @@ FCN_FILE_DIRS += polynomial +plot_PRIVATE_FCN_FILES = \ + polynomial/private/__splinefit__.m + polynomial_FCN_FILES = \ polynomial/compan.m \ polynomial/conv.m \ @@ -24,6 +27,7 @@ polynomial/residue.m \ polynomial/roots.m \ polynomial/spline.m \ + polynomial/splinefit.m \ polynomial/unmkpp.m FCN_FILES += $(polynomial_FCN_FILES)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/private/__splinefit__.m Thu Mar 29 19:13:21 2012 -0400 @@ -0,0 +1,619 @@ +% Copyright (c) 2010, Jonas Lundgren +% All rights reserved. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are +% met: +% +% * Redistributions of source code must retain the above copyright +% notice, this list of conditions and the following disclaimer. +% * Redistributions in binary form must reproduce the above copyright +% notice, this list of conditions and the following disclaimer in +% the documentation and/or other materials provided with the distribution +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. +function pp = __splinefit__(varargin) +%SPLINEFIT Fit a spline to noisy data. +% PP = SPLINEFIT(X,Y,BREAKS) fits a piecewise cubic spline with breaks +% (knots) BREAKS to the noisy data (X,Y). X is a vector and Y is a vector +% or an ND array. If Y is an ND array, then X(j) and Y(:,...,:,j) are +% matched. Use PPVAL to evaluate PP. +% +% PP = SPLINEFIT(X,Y,P) where P is a positive integer interpolates the +% breaks linearly from the sorted locations of X. P is the number of +% spline pieces and P+1 is the number of breaks. +% +% OPTIONAL INPUT +% Argument places 4 to 8 are reserved for optional input. +% These optional arguments can be given in any order: +% +% PP = SPLINEFIT(...,'p') applies periodic boundary conditions to +% the spline. The period length is MAX(BREAKS)-MIN(BREAKS). +% +% PP = SPLINEFIT(...,'r') uses robust fitting to reduce the influence +% from outlying data points. Three iterations of weighted least squares +% are performed. Weights are computed from previous residuals. +% +% PP = SPLINEFIT(...,BETA), where 0 < BETA < 1, sets the robust fitting +% parameter BETA and activates robust fitting ('r' can be omitted). +% Default is BETA = 1/2. BETA close to 0 gives all data equal weighting. +% Increase BETA to reduce the influence from outlying data. BETA close +% to 1 may cause instability or rank deficiency. +% +% PP = SPLINEFIT(...,N) sets the spline order to N. Default is a cubic +% spline with order N = 4. A spline with P pieces has P+N-1 degrees of +% freedom. With periodic boundary conditions the degrees of freedom are +% reduced to P. +% +% PP = SPLINEFIT(...,CON) applies linear constraints to the spline. +% CON is a structure with fields 'xc', 'yc' and 'cc': +% 'xc', x-locations (vector) +% 'yc', y-values (vector or ND array) +% 'cc', coefficients (matrix). +% +% Constraints are linear combinations of derivatives of order 0 to N-2 +% according to +% +% cc(1,j)*y(x) + cc(2,j)*y'(x) + ... = yc(:,...,:,j), x = xc(j). +% +% The maximum number of rows for 'cc' is N-1. If omitted or empty 'cc' +% defaults to a single row of ones. Default for 'yc' is a zero array. +% +% EXAMPLES +% +% % Noisy data +% x = linspace(0,2*pi,100); +% y = sin(x) + 0.1*randn(size(x)); +% % Breaks +% breaks = [0:5,2*pi]; +% +% % Fit a spline of order 5 +% pp = splinefit(x,y,breaks,5); +% +% % Fit a spline of order 3 with periodic boundary conditions +% pp = splinefit(x,y,breaks,3,'p'); +% +% % Constraints: y(0) = 0, y'(0) = 1 and y(3) + y"(3) = 0 +% xc = [0 0 3]; +% yc = [0 1 0]; +% cc = [1 0 1; 0 1 0; 0 0 1]; +% con = struct('xc',xc,'yc',yc,'cc',cc); +% +% % Fit a cubic spline with 8 pieces and constraints +% pp = splinefit(x,y,8,con); +% +% % Fit a spline of order 6 with constraints and periodicity +% pp = splinefit(x,y,breaks,con,6,'p'); +% +% See also SPLINE, PPVAL, PPDIFF, PPINT + +% Author: Jonas Lundgren <splinefit@gmail.com> 2010 + +% 2009-05-06 Original SPLINEFIT. +% 2010-06-23 New version of SPLINEFIT based on B-splines. +% 2010-09-01 Robust fitting scheme added. +% 2010-09-01 Support for data containing NaNs. +% 2011-07-01 Robust fitting parameter added. + +% Check number of arguments +error(nargchk(3,7,nargin)); + +% Check arguments +[x,y,dim,breaks,n,periodic,beta,constr] = arguments(varargin{:}); + +% Evaluate B-splines +base = splinebase(breaks,n); +pieces = base.pieces; +A = ppval(base,x); + +% Bin data +[junk,ibin] = histc(x,[-inf,breaks(2:end-1),inf]); %#ok + +% Sparse system matrix +mx = numel(x); +ii = [ibin; ones(n-1,mx)]; +ii = cumsum(ii,1); +jj = repmat(1:mx,n,1); +if periodic + ii = mod(ii-1,pieces) + 1; + A = sparse(ii,jj,A,pieces,mx); +else + A = sparse(ii,jj,A,pieces+n-1,mx); +end + +% Don't use the sparse solver for small problems +if pieces < 20*n/log(1.7*n) + A = full(A); +end + +% Solve +if isempty(constr) + % Solve Min norm(u*A-y) + u = lsqsolve(A,y,beta); +else + % Evaluate constraints + B = evalcon(base,constr,periodic); + % Solve constraints + [Z,u0] = solvecon(B,constr); + % Solve Min norm(u*A-y), subject to u*B = yc + y = y - u0*A; + A = Z*A; + v = lsqsolve(A,y,beta); + u = u0 + v*Z; +end + +% Periodic expansion of solution +if periodic + jj = mod(0:pieces+n-2,pieces) + 1; + u = u(:,jj); +end + +% Compute polynomial coefficients +ii = [repmat(1:pieces,1,n); ones(n-1,n*pieces)]; +ii = cumsum(ii,1); +jj = repmat(1:n*pieces,n,1); +C = sparse(ii,jj,base.coefs,pieces+n-1,n*pieces); +coefs = u*C; +coefs = reshape(coefs,[],n); + +% Make piecewise polynomial +pp = mkpp(breaks,coefs,dim); + + +%-------------------------------------------------------------------------- +function [x,y,dim,breaks,n,periodic,beta,constr] = arguments(varargin) +%ARGUMENTS Lengthy input checking +% x Noisy data x-locations (1 x mx) +% y Noisy data y-values (prod(dim) x mx) +% dim Leading dimensions of y +% breaks Breaks (1 x (pieces+1)) +% n Spline order +% periodic True if periodic boundary conditions +% beta Robust fitting parameter, no robust fitting if beta = 0 +% constr Constraint structure +% constr.xc x-locations (1 x nx) +% constr.yc y-values (prod(dim) x nx) +% constr.cc Coefficients (?? x nx) + +% Reshape x-data +x = varargin{1}; +mx = numel(x); +x = reshape(x,1,mx); + +% Remove trailing singleton dimensions from y +y = varargin{2}; +dim = size(y); +while numel(dim) > 1 && dim(end) == 1 + dim(end) = []; +end +my = dim(end); + +% Leading dimensions of y +if numel(dim) > 1 + dim(end) = []; +else + dim = 1; +end + +% Reshape y-data +pdim = prod(dim); +y = reshape(y,pdim,my); + +% Check data size +if mx ~= my + mess = 'Last dimension of array y must equal length of vector x.'; + error('arguments:datasize',mess) +end + +% Treat NaNs in x-data +inan = find(isnan(x)); +if ~isempty(inan) + x(inan) = []; + y(:,inan) = []; + mess = 'All data points with NaN as x-location will be ignored.'; + warning('arguments:nanx',mess) +end + +% Treat NaNs in y-data +inan = find(any(isnan(y),1)); +if ~isempty(inan) + x(inan) = []; + y(:,inan) = []; + mess = 'All data points with NaN in their y-value will be ignored.'; + warning('arguments:nany',mess) +end + +% Check number of data points +mx = numel(x); +if mx == 0 + error('arguments:nodata','There must be at least one data point.') +end + +% Sort data +if any(diff(x) < 0) + [x,isort] = sort(x); + y = y(:,isort); +end + +% Breaks +if isscalar(varargin{3}) + % Number of pieces + p = varargin{3}; + if ~isreal(p) || ~isfinite(p) || p < 1 || fix(p) < p + mess = 'Argument #3 must be a vector or a positive integer.'; + error('arguments:breaks1',mess) + end + if x(1) < x(end) + % Interpolate breaks linearly from x-data + dx = diff(x); + ibreaks = linspace(1,mx,p+1); + [junk,ibin] = histc(ibreaks,[0,2:mx-1,mx+1]); %#ok + breaks = x(ibin) + dx(ibin).*(ibreaks-ibin); + else + breaks = x(1) + linspace(0,1,p+1); + end +else + % Vector of breaks + breaks = reshape(varargin{3},1,[]); + if isempty(breaks) || min(breaks) == max(breaks) + mess = 'At least two unique breaks are required.'; + error('arguments:breaks2',mess); + end +end + +% Unique breaks +if any(diff(breaks) <= 0) + breaks = unique(breaks); +end + +% Optional input defaults +n = 4; % Cubic splines +periodic = false; % No periodic boundaries +robust = false; % No robust fitting scheme +beta = 0.5; % Robust fitting parameter +constr = []; % No constraints + +% Loop over optional arguments +for k = 4:nargin + a = varargin{k}; + if ischar(a) && isscalar(a) && lower(a) == 'p' + % Periodic conditions + periodic = true; + elseif ischar(a) && isscalar(a) && lower(a) == 'r' + % Robust fitting scheme + robust = true; + elseif isreal(a) && isscalar(a) && isfinite(a) && a > 0 && a < 1 + % Robust fitting parameter + beta = a; + robust = true; + elseif isreal(a) && isscalar(a) && isfinite(a) && a > 0 && fix(a) == a + % Spline order + n = a; + elseif isstruct(a) && isscalar(a) + % Constraint structure + constr = a; + else + error('arguments:nonsense','Failed to interpret argument #%d.',k) + end +end + +% No robust fitting +if ~robust + beta = 0; +end + +% Check exterior data +h = diff(breaks); +xlim1 = breaks(1) - 0.01*h(1); +xlim2 = breaks(end) + 0.01*h(end); +if x(1) < xlim1 || x(end) > xlim2 + if periodic + % Move data inside domain + P = breaks(end) - breaks(1); + x = mod(x-breaks(1),P) + breaks(1); + % Sort + [x,isort] = sort(x); + y = y(:,isort); + else + mess = 'Some data points are outside the spline domain.'; + warning('arguments:exteriordata',mess) + end +end + +% Return +if isempty(constr) + return +end + +% Unpack constraints +xc = []; +yc = []; +cc = []; +names = fieldnames(constr); +for k = 1:numel(names) + switch names{k} + case {'xc'} + xc = constr.xc; + case {'yc'} + yc = constr.yc; + case {'cc'} + cc = constr.cc; + otherwise + mess = 'Unknown field ''%s'' in constraint structure.'; + warning('arguments:unknownfield',mess,names{k}) + end +end + +% Check xc +if isempty(xc) + mess = 'Constraints contains no x-locations.'; + error('arguments:emptyxc',mess) +else + nx = numel(xc); + xc = reshape(xc,1,nx); +end + +% Check yc +if isempty(yc) + % Zero array + yc = zeros(pdim,nx); +elseif numel(yc) == 1 + % Constant array + yc = zeros(pdim,nx) + yc; +elseif numel(yc) ~= pdim*nx + % Malformed array + error('arguments:ycsize','Cannot reshape yc to size %dx%d.',pdim,nx) +else + % Reshape array + yc = reshape(yc,pdim,nx); +end + +% Check cc +if isempty(cc) + cc = ones(size(xc)); +elseif numel(size(cc)) ~= 2 + error('arguments:ccsize1','Constraint coefficients cc must be 2D.') +elseif size(cc,2) ~= nx + mess = 'Last dimension of cc must equal length of xc.'; + error('arguments:ccsize2',mess) +end + +% Check high order derivatives +if size(cc,1) >= n + if any(any(cc(n:end,:))) + mess = 'Constraints involve derivatives of order %d or larger.'; + error('arguments:difforder',mess,n-1) + end + cc = cc(1:n-1,:); +end + +% Check exterior constraints +if min(xc) < xlim1 || max(xc) > xlim2 + if periodic + % Move constraints inside domain + P = breaks(end) - breaks(1); + xc = mod(xc-breaks(1),P) + breaks(1); + else + mess = 'Some constraints are outside the spline domain.'; + warning('arguments:exteriorconstr',mess) + end +end + +% Pack constraints +constr = struct('xc',xc,'yc',yc,'cc',cc); + + +%-------------------------------------------------------------------------- +function pp = splinebase(breaks,n) +%SPLINEBASE Generate B-spline base PP of order N for breaks BREAKS + +breaks = breaks(:); % Breaks +breaks0 = breaks'; % Initial breaks +h = diff(breaks); % Spacing +pieces = numel(h); % Number of pieces +deg = n - 1; % Polynomial degree + +% Extend breaks periodically +if deg > 0 + if deg <= pieces + hcopy = h; + else + hcopy = repmat(h,ceil(deg/pieces),1); + end + % to the left + hl = hcopy(end:-1:end-deg+1); + bl = breaks(1) - cumsum(hl); + % and to the right + hr = hcopy(1:deg); + br = breaks(end) + cumsum(hr); + % Add breaks + breaks = [bl(deg:-1:1); breaks; br]; + h = diff(breaks); + pieces = numel(h); +end + +% Initiate polynomial coefficients +coefs = zeros(n*pieces,n); +coefs(1:n:end,1) = 1; + +% Expand h +ii = [1:pieces; ones(deg,pieces)]; +ii = cumsum(ii,1); +ii = min(ii,pieces); +H = h(ii(:)); + +% Recursive generation of B-splines +for k = 2:n + % Antiderivatives of splines + for j = 1:k-1 + coefs(:,j) = coefs(:,j).*H/(k-j); + end + Q = sum(coefs,2); + Q = reshape(Q,n,pieces); + Q = cumsum(Q,1); + c0 = [zeros(1,pieces); Q(1:deg,:)]; + coefs(:,k) = c0(:); + % Normalize antiderivatives by max value + fmax = repmat(Q(n,:),n,1); + fmax = fmax(:); + for j = 1:k + coefs(:,j) = coefs(:,j)./fmax; + end + % Diff of adjacent antiderivatives + coefs(1:end-deg,1:k) = coefs(1:end-deg,1:k) - coefs(n:end,1:k); + coefs(1:n:end,k) = 0; +end + +% Scale coefficients +scale = ones(size(H)); +for k = 1:n-1 + scale = scale./H; + coefs(:,n-k) = scale.*coefs(:,n-k); +end + +% Reduce number of pieces +pieces = pieces - 2*deg; + +% Sort coefficients by interval number +ii = [n*(1:pieces); deg*ones(deg,pieces)]; +ii = cumsum(ii,1); +coefs = coefs(ii(:),:); + +% Make piecewise polynomial +pp = mkpp(breaks0,coefs,n); + + +%-------------------------------------------------------------------------- +function B = evalcon(base,constr,periodic) +%EVALCON Evaluate linear constraints + +% Unpack structures +breaks = base.breaks; +pieces = base.pieces; +n = base.order; +xc = constr.xc; +cc = constr.cc; + +% Bin data +[junk,ibin] = histc(xc,[-inf,breaks(2:end-1),inf]); %#ok + +% Evaluate constraints +nx = numel(xc); +B0 = zeros(n,nx); +for k = 1:size(cc,1) + if any(cc(k,:)) + B0 = B0 + repmat(cc(k,:),n,1).*ppval(base,xc); + end + % Differentiate base + coefs = base.coefs(:,1:n-k); + for j = 1:n-k-1 + coefs(:,j) = (n-k-j+1)*coefs(:,j); + end + base.coefs = coefs; + base.order = n-k; +end + +% Sparse output +ii = [ibin; ones(n-1,nx)]; +ii = cumsum(ii,1); +jj = repmat(1:nx,n,1); +if periodic + ii = mod(ii-1,pieces) + 1; + B = sparse(ii,jj,B0,pieces,nx); +else + B = sparse(ii,jj,B0,pieces+n-1,nx); +end + + +%-------------------------------------------------------------------------- +function [Z,u0] = solvecon(B,constr) +%SOLVECON Find a particular solution u0 and null space Z (Z*B = 0) +% for constraint equation u*B = yc. + +yc = constr.yc; +tol = 1000*eps; + +% Remove blank rows +ii = any(B,2); +B2 = full(B(ii,:)); + +% Null space of B2 +if isempty(B2) + Z2 = []; +else + % QR decomposition with column permutation + [Q,R,dummy] = qr(B2); %#ok + R = abs(R); + jj = all(R < R(1)*tol, 2); + Z2 = Q(:,jj)'; +end + +% Sizes +[m,ncon] = size(B); +m2 = size(B2,1); +nz = size(Z2,1); + +% Sparse null space of B +Z = sparse(nz+1:nz+m-m2,find(~ii),1,nz+m-m2,m); +Z(1:nz,ii) = Z2; + +% Warning rank deficient +if nz + ncon > m2 + mess = 'Rank deficient constraints, rank = %d.'; + warning('solvecon:deficient',mess,m2-nz); +end + +% Particular solution +u0 = zeros(size(yc,1),m); +if any(yc(:)) + % Non-homogeneous case + u0(:,ii) = yc/B2; + % Check solution + if norm(u0*B - yc,'fro') > norm(yc,'fro')*tol + mess = 'Inconsistent constraints. No solution within tolerance.'; + error('solvecon:inconsistent',mess) + end +end + + +%-------------------------------------------------------------------------- +function u = lsqsolve(A,y,beta) +%LSQSOLVE Solve Min norm(u*A-y) + +% Avoid sparse-complex limitations +if issparse(A) && ~isreal(y) + A = full(A); +end + +% Solution +u = y/A; + +% Robust fitting +if beta > 0 + [m,n] = size(y); + alpha = 0.5*beta/(1-beta)/m; + for k = 1:3 + % Residual + r = u*A - y; + rr = r.*conj(r); + rrmean = sum(rr,2)/n; + rrmean(~rrmean) = 1; + rrhat = (alpha./rrmean)'*rr; + % Weights + w = exp(-rrhat); + spw = spdiags(w',0,n,n); + % Solve weighted problem + u = (y*spw)/(A*spw); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/splinefit.m Thu Mar 29 19:13:21 2012 -0400 @@ -0,0 +1,238 @@ +## Copyright (C) 2012 Ben Abbott, Jonas Lundgren +## +## 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 3 of the License, 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, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{pp} =} splinefit (@var{x}, @var{y}, @var{breaks}) +## Fits a piecewise cubic spline with breaks (knots) @var{breaks} to the +## noisy data, @var{x} and @var{y}. @var{x} is a vector, and @var{y} +## a vector or ND array. If @var{y} is an ND array, then @var{x}(j) +## is matched to @var{y}(:,...,:,j). +## +## The fitted spline is returned as a piece-wise polynomial, @var{pp}, and +## may be evaluated using @code{ppval}. +## +## @deftypefnx {Function File} {@var{pp} =} splinefit (@var{x}, @var{y}, @var{p}) +## @var{p} is a positive integer defining the number of linearly spaced +## intervals along @var{x}. @var{p} is the number of intervalas and +## @var{p}+1 is the number of breaks. +## +## @deftypefnx {Function File} {@var{pp} =} splinefit (@dots{}, "periodic", @var{periodic}) +## @deftypefnx {Function File} {@var{pp} =} splinefit (@dots{}, "robust", @var{robust}) +## @deftypefnx {Function File} {@var{pp} =} splinefit (@dots{}, "beta", @var{beta}) +## @deftypefnx {Function File} {@var{pp} =} splinefit (@dots{}, "order", @var{order}) +## @deftypefnx {Function File} {@var{pp} =} splinefit (@dots{}, "constraints", @var{constraints}) +## +## The optional property @var{periodic} is a logical value which specifies +## whether a periodic boundary condition is applied to the spline. The +## length of the period is @code{max(@var{breaks})-min(@var{breaks})}. +## The default value is @code{false}. +## +## The optional property @var{robust} is a logical value which specifies +## if robust fitting is to be applied to reduce the influence of outlying +## data points. Three iterations of weighted least squares are performed. +## Weights are computed from previous residuals. The sensitivity of outlier +## identification is controlled by the property @var{beta}. The value of +## @var{beta} is stricted to the range, 0 < @var{beta} < 1. The default +## value is @var{beta} = 1/2. Values close to 0 give all data equal +## weighting. Increasing values of @var{beta} reduce the influence of +## outlying data. Values close to unity may cause instability or rank +## deficiency. +## +## @var{order} sets the order polynomials used to construct the spline. +## The default is a cubic, @var{order}=3. A spline with P pieces has +## P+@var{order} degrees of freedom. With periodic boundary conditions +## the degrees of freedom are reduced to P. +## +## The optional property, @var{constaints}, is a structure specifying +## linear constraints on the fit. The structure has three fields, "xc", +## "yc", and "cc". +## +## @table @asis +## @item "xc" +## x-locations of the constraints (vector) with a size identical to @var{x}. +## @item "yc" +## Constaining values with a size identical to @var{y}. The default +## is an array of zeros. +## @item "cc" +## Coefficients (matrix). The default is an array of ones. The number of +## rows is limited to the order of the piece-wise polynomials, @var{order}. +## @end table +## +## Constraints are linear combinations of derivatives of order 0 to +## @var{order}-1 according to +## +## @example +## @group +## @tex +## $cc(1,j) \cdot y(x) + cc(2,j) \cdot y\prime(x) + ... = yc(:,\dots,:,j), \quad x = xc(j)$. +## @end tex +## @ifnottex +## cc(1,j) * y(x) + cc(2,j) * y'(x) + ... = yc(:,...,:,j), x = xc(j). +## @end ifnottex +## @end group +## @end example +## +## @seealso{interp1, unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps} +## @end deftypefn + +%!demo +%! % Noisy data +%! x = linspace (0, 2*pi, 100); +%! y = sin (x) + 0.1 * randn (size (x)); +%! % Breaks +%! breaks = [0:5, 2*pi]; +%! % Fit a spline of order 5 +%! pp = splinefit (x, y, breaks, "order", 4); +%! clf () +%! plot (x, y, "s", x, ppval (pp, x), "r", breaks, ppval (pp, breaks), "+r") +%! xlabel ("Independent Variable") +%! ylabel ("Dependent Variable") +%! title ("Fit a piece-wise polynomial of order 4"); +%! legend ({"data", "fit", "breaks"}) +%! axis tight +%! ylim auto + +%!demo +%! % Noisy data +%! x = linspace (0,2*pi, 100); +%! y = sin (x) + 0.1 * randn (size (x)); +%! % Breaks +%! breaks = [0:5, 2*pi]; +%! % Fit a spline of order 3 with periodic boundary conditions +%! pp = splinefit (x, y, breaks, "order", 2, "periodic", true); +%! clf () +%! plot (x, y, "s", x, ppval (pp, x), "r", breaks, ppval (pp, breaks), "+r") +%! xlabel ("Independent Variable") +%! ylabel ("Dependent Variable") +%! title ("Fit a periodic piece-wise polynomial of order 2"); +%! legend ({"data", "fit", "breaks"}) +%! axis tight +%! ylim auto + +%!demo +%! % Noisy data +%! x = linspace (0, 2*pi, 100); +%! y = sin (x) + 0.1 * randn (size (x)); +%! % Breaks +%! breaks = [0:5, 2*pi]; +%! % Constraints: y(0) = 0, y"(0) = 1 and y(3) + y"(3) = 0 +%! xc = [0 0 3]; +%! yc = [0 1 0]; +%! cc = [1 0 1; 0 1 0; 0 0 1]; +%! con = struct ("xc", xc, "yc", yc, "cc", cc); +%! % Fit a cubic spline with 8 pieces and constraints +%! pp = splinefit (x, y, 8, "constraints", con); +%! clf () +%! plot (x, y, "s", x, ppval (pp, x), "r", breaks, ppval (pp, breaks), "+r") +%! xlabel ("Independent Variable") +%! ylabel ("Dependent Variable") +%! title ("Fit a cubic spline with constraints") +%! legend ({"data", "fit", "breaks"}) +%! axis tight +%! ylim auto + +%!demo +%! % Noisy data +%! x = linspace (0, 2*pi, 100); +%! y = sin (x) + 0.1 * randn (size (x)); +%! % Breaks +%! breaks = [0:5, 2*pi]; +%! xc = [0 0 3]; +%! yc = [0 1 0]; +%! cc = [1 0 1; 0 1 0; 0 0 1]; +%! con = struct ("xc", xc, "yc", yc, "cc", cc); +%! % Fit a spline of order 6 with constraints and periodicity +%! pp = splinefit (x, y, breaks, "constraints", con, "order", 5, "periodic", true); +%! clf () +%! plot (x, y, "s", x, ppval (pp, x), "r", breaks, ppval (pp, breaks), "+r") +%! xlabel ("Independent Variable") +%! ylabel ("Dependent Variable") +%! title ("Fit a 5th order piece-wise periodic polynomial with constraints") +%! legend ({"data", "fit", "breaks"}) +%! axis tight +%! ylim auto + +function pp = splinefit (x, y, breaks, varargin) + if (nargin > 3) + n = cellfun (@ischar, varargin, "uniformoutput", true); + varargin(n) = lower (varargin(n)); + try + props = struct (varargin{:}); + catch + print_usage (); + end_try_catch + else + props = struct (); + endif + fields = fieldnames (props); + for f = 1:numel(fields) + if (! any (strcmp (fields{f}, {"periodic", "robust", "beta", ... + "order", "constraints"}))) + error (sprintf ("%s:invalidproperty", mfilename ()), + sprintf ("""%s"" is not recongizied", fields{f})) + endif + endfor + args = {}; + if (isfield (props, "periodic") && props.periodic) + args{end+1} = "p"; + endif + if (isfield (props, "robust") && props.robust) + args{end+1} = "r"; + endif + if (isfield (props, "beta")) + if (0 < props.beta && props.beta < 1) + args{end+1} = props.beta; + else + error (sprintf ("%s:invalidbeta", mfilename), + "Invalid beta parameter (0 < beta < 1)") + endif + endif + if (isfield (props, "order")) + if (props.order >= 0) + args{end+1} = props.order + 1; + else + error (sprintf ("%s:invalidorder", mfilename), "Invalid order") + endif + endif + if (isfield (props, "constraints")) + args{end+1} = props.constraints; + endif + if (nargin < 3) + print_usage (); + elseif (! isnumeric (breaks) || ! isvector (breaks)) + print_usage (); + endif + pp = __splinefit__ (x, y, breaks, args{:}); +endfunction + +%!shared xb, yb, x +%! xb = 0:2:10; +%! yb = randn (size (xb)); +%! x = 0:0.1:10; + +%!test +%! y = interp1 (xb, yb, x, "linear"); +%! assert (ppval (splinefit (x, y, xb, "order", 1), x), y, 10 * eps ()); +%!test +%! y = interp1 (xb, yb, x, "spline"); +%! assert (ppval (splinefit (x, y, xb, "order", 3), x), y, 10 * eps ()); +%!test +%! y = interp1 (xb, yb, x, "spline"); +%! assert (ppval (splinefit (x, y, xb), x), y, 10 * eps ()); + +