1 % Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
3 % Author: Jaroslav Hajek <highegg@gmail.com>
5 % This file is part of OctGPR.
7 % OctGPR is free software; you can redistribute it and/or modify
8 % it under the terms of the GNU General Public License as published by
9 % the Free Software Foundation; either version 3 of the License, or
10 % (at your option) any later version.
12 % This program is distributed in the hope that it will be useful,
13 % but WITHOUT ANY WARRANTY; without even the implied warranty of
14 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 % GNU General Public License for more details.
17 % You should have received a copy of the GNU General Public License
18 % along with this software; see the file COPYING. If not, see
19 % <http://www.gnu.org/licenses/>.
23 % @deftypefn {Function File} demo_octgpr (1, nsamp = 150)
24 % @deftypefnx {Function File} demo_octgpr (2, ncnt = 20, npt = 500)
25 % @deftypefnx {Function File} demo_octgpr (3, ncnt = 50, nsamp = 500)
26 % OctGPR package demo function.
27 % First argument selects available demos:
30 % @item 1. GPR regression demo @*
31 % A function is sampled (with small noise), then reconstructed using GPR
32 % regression. @var{nsamp} specifies the number of samples.
33 % @seealso{gpr_train, gpr_predict}
34 % @item 2. RBF centers selection demo @*
35 % Radial basis centers are selected amongst random points.
36 % @var{ncnt} specifies number of centers, @var{npt} number of points.
37 % @item 2. PGP regression demo @*
38 % A function is densely sampled (with small noise),
39 % radial basis centers are selected, then the function is reconstructed
40 % using PGP regression. @var{nsamp} specifies the number of samples,
41 % @var{ncnt} specifies number of centers.
44 function demo_octgpr (number, varargin)
50 elseif (ischar (number))
52 elseif (isscalar (number))
54 if (! isempty (prntfmt))
55 figure (gcf, "visible", "off");
60 demo_octgpr1 (varargin{:})
62 demo_octgpr2 (varargin{:})
64 demo_octgpr3 (varargin{:})
66 error ("demo_octgpr: invalid demo number")
74 function demo_octgpr_pause (idemo, iplot)
76 if (isempty (prntfmt))
79 print (sprintf (prntfmt, idemo, iplot));
84 % define the test function (the well-known matlab "peaks" plus some sines)
85 function z = testfun1 (x, y)
86 z = 4 + 3 * (1-x).^2 .* exp(-(x.^2) - (y+1).^2) ...
87 + 10 * (x/5 - x.^3 - y.^5) .* exp(-x.^2 - y.^2) ...
88 - 1/3 * exp(-(x+1).^2 - y.^2) ...
89 + 2*sin (x + y + 1e-1*x.*y);
92 function demo_octgpr1 (nsamp = 150)
94 tit = "a peaked surface";
97 % create the mesh onto which to interpolate
98 t = linspace (-3, 3, 50);
99 [xi,yi] = meshgrid (t, t);
102 zi = testfun1 (xi, yi);
103 zimax = max (vec (zi)); zimin = min (vec (zi));
108 contourf (xi, yi, zi, 20);
109 demo_octgpr_pause (1, 1);
111 tit = sprintf ("sampled at %d random points", nsamp);
113 % create random samples
114 xs = rand (nsamp,1); ys = rand (nsamp,1);
115 xs = 6*xs-3; ys = 6*ys - 3;
116 % evaluate at random samples
117 zs = testfun1 (xs, ys);
121 plot3 (xs, ys, zs, ".+");
129 demo_octgpr_pause (1, 2);
131 tit = "GPR model with heuristic hypers";
133 ths = 1 ./ std (xys);
134 GPM = gpr_train (xys, zs, ths, 1e-5);
135 zm = gpr_predict (GPM, [vec(xi) vec(yi)]);
136 zm = reshape (zm, size(zi));
137 zm = min (zm, zimax); zm = max (zm, zimin);
143 contourf (xi, yi, zm, 20);
146 demo_octgpr_pause (1, 3);
148 tit = "GPR model with MLE training";
151 GPM = gpr_train (xys, zs, ths, 1e-5, {"tol", 1e-5, "maxev", 400, "numin", 1e-8});
152 zm = gpr_predict (GPM, [vec(xi) vec(yi)]);
153 zm = reshape (zm, size (zi));
154 zm = min (zm, zimax); zm = max (zm, zimin);
160 contourf (xi, yi, zm, 20);
163 demo_octgpr_pause (1, 4);
168 function demo_octgpr2 (ncnt = 50, npt = 500)
171 npt = ncnt*ceil (npt/ncnt);
173 cs = min (pdist2_mw (U, 2) + diag (Inf (ncnt, 1)));
174 X = repmat (U, npt/ncnt, 1) + repmat (cs', npt/ncnt, 2) .* randn (npt, 2);
175 disp ("slightly clustered random points")
176 plot (X(:,1), X(:,2), "+");
177 demo_octgpr_pause (2, 1);
179 [U, ur] = rbf_centers(X, ncnt);
181 fi = linspace (0, 2*pi, 20);
182 ncolors = rows (colormap);
185 xc = U(i,1) + ur(i) * cos (fi);
186 yc = U(i,2) + ur(i) * sin (fi);
190 demo_octgpr_pause (2, 2);
195 function demo_octgpr3 (ncnt = 100, nsamp = 1000)
198 tit = "a peaked surface";
201 % create the mesh onto which to interpolate
202 t = linspace (-3, 3, 50);
203 [xi,yi] = meshgrid (t, t);
206 zi = testfun1 (xi, yi);
207 zimax = max (vec (zi)); zimin = min (vec (zi));
212 contourf (xi, yi, zi, 20);
213 demo_octgpr_pause (1, 1);
215 tit = sprintf ("sampled at %d random points, selected %d centers", nsamp, ncnt);
217 % create random samples
218 xs = rand (nsamp,1); ys = rand (nsamp,1);
219 xs = 6*xs-3; ys = 6*ys - 3;
220 % evaluate at random samples
221 zs = testfun1 (xs, ys);
224 % select centers using k-means
225 xyc = rbf_centers (xys, ncnt);
226 xc = xyc(:,1); yc = xyc(:,2);
229 plot3 (xs, ys, zs, ".+");
240 demo_octgpr_pause (1, 2);
242 tit = "PGP model with heuristic hypers";
244 ths = 1 ./ std (xyc);
245 GPM = pgp_train (xys, xyc, zs, ths, 1e-5);
246 zm = pgp_predict (GPM, [vec(xi) vec(yi)]);
247 zm = reshape (zm, size(zi));
248 zm = min (zm, zimax); zm = max (zm, zimin);
254 contourf (xi, yi, zm, 20);
258 demo_octgpr_pause (1, 3);
260 tit = "PGP model with MLE training";
263 GPM = pgp_train (xys, xyc, zs, ths, 1e-3, {"tol", 1e-5, "maxev", 400});
264 zm = pgp_predict (GPM, [vec(xi) vec(yi)]);
265 zm = reshape (zm, size (zi));
266 zm = min (zm, zimax); zm = max (zm, zimin);
272 contourf (xi, yi, zm, 20);
276 demo_octgpr_pause (1, 4);