|
14509
|
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 |
