Mercurial > hg > octave-nkf

comparison doc/interpreter/splineimages.m @ 14509:a88f8e4fae56

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
children e12945668746
comparison
equal deleted inserted replaced
14508:0901f926ed50 14509:a88f8e4fae56
1 ## Copyright (C) 2012 Ben Abbott, Jonas Lundgren
2 ##
3 ## This file is part of Octave.
4 ##
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
9 ##
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
18
19 function splineimages (nm, typ)
20 graphics_toolkit ("gnuplot");
21 set_print_size ();
22 hide_output ();
23 if (strcmp (typ, "png"))
24 set (0, "defaulttextfontname", "*");
25 endif
26 if (strcmp (typ, "eps"))
27 d_typ = "-depsc2";
28 else
29 d_typ = cstrcat ("-d", typ);
30 endif
31
32 if (strcmp (typ, "txt"))
33 image_as_txt (nm);
34 elseif (strcmp (nm, "splinefit1")) ## Breaks and Pieces
35 x = 2 * pi * rand (1, 200);
36 y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
37 ## Uniform breaks
38 breaks = linspace (0, 2 * pi, 41); ## 41 breaks, 40 pieces
39 pp1 = splinefit (x, y, breaks);
40 ## Breaks interpolated from data
41 pp2 = splinefit (x, y, 10); ## 11 breaks, 10 pieces
42 ## Plot
43 xx = linspace (0, 2 * pi, 400);
44 y1 = ppval (pp1, xx);
45 y2 = ppval (pp2, xx);
46 plot (x, y, ".", xx, [y1; y2])
47 axis tight
48 ylim ([-2.5 2.5])
49 legend ("data", "41 breaks, 40 pieces", "11 breaks, 10 pieces")
50 print (cstrcat (nm, ".", typ), d_typ)
51 elseif (strcmp (nm, "splinefit2")) ## Spline orders
52 ## Data (200 points)
53 x = 2 * pi * rand (1, 200);
54 y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
55 ## Splines
56 pp1 = splinefit (x, y, 8, "order", 0); ## Piecewise constant
57 pp2 = splinefit (x, y, 8, "order", 1); ## Piecewise linear
58 pp3 = splinefit (x, y, 8, "order", 2); ## Piecewise quadratic
59 pp4 = splinefit (x, y, 8, "order", 3); ## Piecewise cubic
60 pp5 = splinefit (x, y, 8, "order", 4); ## Etc.
61 ## Plot
62 xx = linspace (0, 2 * pi, 400);
63 y1 = ppval (pp1, xx);
64 y2 = ppval (pp2, xx);
65 y3 = ppval (pp3, xx);
66 y4 = ppval (pp4, xx);
67 y5 = ppval (pp5, xx);
68 plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
69 axis tight
70 ylim ([-2.5 2.5])
71 legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})
72 print (cstrcat (nm, ".", typ), d_typ)
73 elseif (strcmp (nm, "splinefit3"))
74 ## Data (100 points)
75 x = 2 * pi * [0, (rand (1, 98)), 1];
76 y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
77 ## No constraints
78 pp1 = splinefit (x, y, 10, "order", 5);
79 ## Periodic boundaries
80 pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
81 ## Plot
82 xx = linspace (0, 2 * pi, 400);
83 y1 = ppval (pp1, xx);
84 y2 = ppval (pp2, xx);
85 plot (x, y, ".", xx, [y1; y2])
86 axis tight
87 ylim ([-2 3])
88 legend ({"data", "no constraints", "periodic"})
89 print (cstrcat (nm, ".", typ), d_typ)
90 elseif (strcmp (nm, "splinefit4"))
91 ## Data (200 points)
92 x = 2 * pi * rand (1, 200);
93 y = sin (2 * x) + 0.1 * randn (size (x));
94 ## Breaks
95 breaks = linspace (0, 2 * pi, 10);
96 ## Clamped endpoints, y = y" = 0
97 xc = [0, 0, 2*pi, 2*pi];
98 cc = [(eye (2)), (eye (2))];
99 con = struct ("xc", xc, "cc", cc);
100 pp1 = splinefit (x, y, breaks, "constraints", con);
101 ## Hinged periodic endpoints, y = 0
102 con = struct ("xc", 0);
103 pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
104 ## Plot
105 xx = linspace (0, 2 * pi, 400);
106 y1 = ppval (pp1, xx);
107 y2 = ppval (pp2, xx);
108 plot (x, y, ".", xx, [y1; y2])
109 axis tight
110 ylim ([-1.5 1.5])
111 legend({"data", "clamped", "hinged periodic"})
112 print (cstrcat (nm, ".", typ), d_typ)
113 elseif (strcmp (nm, "splinefit5"))
114 ## Truncated data
115 x = [0, 1, 2, 4, 8, 16, 24, 40, 56, 72, 80] / 80;
116 y = [0, 28, 39, 53, 70, 86, 90, 79, 55, 22, 2] / 1000;
117 xy = [x; y];
118 ## Curve length parameter
119 ds = sqrt (diff (x).^2 + diff (y).^2);
120 s = [0, cumsum(ds)];
121 ## Constraints at s = 0: (x,y) = (0,0), (dx/ds,dy/ds) = (0,1)
122 con = struct ("xc", [0 0], "yc", [0 0; 0 1], "cc", eye (2));
123 ## Fit a spline with 4 pieces
124 pp = splinefit (s, xy, 4, "constraints", con);
125 ## Plot
126 ss = linspace (0, s(end), 400);
127 xyfit = ppval (pp, ss);
128 xyb = ppval(pp, pp.breaks);
129 plot (x, y, ".", xyfit(1,:), xyfit(2,:), "r", xyb(1,:), xyb(2,:), "ro")
130 legend ({"data", "spline", "breaks"})
131 axis tight
132 ylim ([0 0.1])
133 print (cstrcat (nm, ".", typ), d_typ)
134 elseif (strcmp (nm, "splinefit6"))
135 ## Data
136 x = linspace (0, 2*pi, 200);
137 y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
138 ## Add outliers
139 x = [x, linspace(0,2*pi,60)];
140 y = [y, -ones(1,60)];
141 ## Fit splines with hinged conditions
142 con = struct ("xc", [0, 2*pi]);
143 pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25); ## Robust fitting
144 pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75); ## Robust fitting
145 pp3 = splinefit (x, y, 8, "constraints", con); ## No robust fitting
146 ## Plot
147 xx = linspace (0, 2*pi, 400);
148 y1 = ppval (pp1, xx);
149 y2 = ppval (pp2, xx);
150 y3 = ppval (pp3, xx);
151 plot (x, y, ".", xx, [y1; y2; y3])
152 legend({"data with outliers","robust, beta = 0.25", ...
153 "robust, beta = 0.75", "no robust fitting"})
154 axis tight
155 ylim ([-2 2])
156 print (cstrcat (nm, ".", typ), d_typ)
157 endif
158 hide_output ();
159 endfunction
160
161 function set_print_size ()
162 image_size = [5.0, 3.5]; # in inches, 16:9 format
163 border = 0; # For postscript use 50/72
164 set (0, "defaultfigurepapertype", "<custom>");
165 set (0, "defaultfigurepaperorientation", "landscape");
166 set (0, "defaultfigurepapersize", image_size + 2*border);
167 set (0, "defaultfigurepaperposition", [border, border, image_size]);
168 endfunction
169
170 ## Use this function before plotting commands and after every call to
171 ## print since print() resets output to stdout (unfortunately, gnpulot
172 ## can't pop output as it can the terminal type).
173 function hide_output ()
174 f = figure (1);
175 set (f, "visible", "off");
176 endfunction
177
178 ## generate something for the texinfo @image command to process
179 function image_as_txt(nm)
180 fid = fopen (sprintf ("##s.txt", nm), "wt");
181 fputs (fid, "\n");
182 fputs (fid, "+---------------------------------+\n");
183 fputs (fid, "| Image unavailable in text mode. |\n");
184 fputs (fid, "+---------------------------------+\n");
185 fclose (fid);
186 endfunction
187
188 %!demo
189 %! for s = 1:6
190 %! splineimages (sprintf ("splinefit##d", s), "pdf")
191 %! endfor
192