octave_packages/octgpr-1.2.0/demo_octgpr.m

   1 % Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
   2 %
   3 % Author: Jaroslav Hajek <highegg@gmail.com>
   4 %
   5 % This file is part of OctGPR.
   6 %
   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.
  11 %
  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.
  16 %
  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/>.
  20 %
  21
  22 % -*- texinfo -*-
  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:
  28 %
  29 % @itemize
  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.
  42 % @end itemize
  43 % @end deftypefn
  44 function demo_octgpr (number, varargin)
  45
  46   global prntfmt = '';
  47
  48   if (nargin < 1)
  49     print_usage ();
  50   elseif (ischar (number))
  51     prntfmt = number;
  52   elseif (isscalar (number))
  53     figure ();
  54     if (! isempty (prntfmt))
  55       figure (gcf, "visible", "off");
  56     endif
  57
  58     switch (number)
  59     case 1
  60       demo_octgpr1 (varargin{:})
  61     case 2
  62       demo_octgpr2 (varargin{:})
  63     case 3
  64       demo_octgpr3 (varargin{:})
  65     otherwise
  66       error ("demo_octgpr: invalid demo number")
  67     endswitch
  68   else
  69     print_usage ();
  70   endif
  71
  72 endfunction
  73
  74 function demo_octgpr_pause (idemo, iplot)
  75   global prntfmt;
  76   if (isempty (prntfmt))
  77     pause;
  78   else
  79     print (sprintf (prntfmt, idemo, iplot));
  80     fflush (stdout);
  81   endif
  82 endfunction
  83
  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);
  90 endfunction
  91
  92 function demo_octgpr1 (nsamp = 150)
  93   global prntfmt;
  94   tit = "a peaked surface";
  95   disp (tit);
  96
  97   % create the mesh onto which to interpolate
  98   t = linspace (-3, 3, 50);
  99   [xi,yi] = meshgrid (t, t);
 100
 101   % evaluate
 102   zi = testfun1 (xi, yi);
 103   zimax = max (vec (zi)); zimin = min (vec (zi));
 104   subplot (2, 2, 1);
 105   mesh (xi, yi, zi);
 106   title (tit);
 107   subplot (2, 2, 3);
 108   contourf (xi, yi, zi, 20);
 109   demo_octgpr_pause (1, 1);
 110
 111   tit = sprintf ("sampled at %d random points", nsamp);
 112   disp (tit);
 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);
 118   xys = [xs ys];
 119
 120   subplot (2, 2, 2);
 121   plot3 (xs, ys, zs, ".+");
 122   title (tit);
 123   subplot (2, 2, 3);
 124   hold on
 125   plot (xs, ys, "+6");
 126   hold off
 127   subplot (2, 2, 4);
 128   plot (xs, ys, ".+");
 129   demo_octgpr_pause (1, 2);
 130
 131   tit = "GPR model with heuristic hypers";
 132   disp (tit);
 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);
 138   subplot (2, 2, 2);
 139   mesh (xi, yi, zm);
 140   title (tit);
 141   subplot(2, 2, 4);
 142   hold on
 143   contourf (xi, yi, zm, 20);
 144   plot (xs, ys, "+6");
 145   hold off
 146   demo_octgpr_pause (1, 3);
 147
 148   tit = "GPR model with MLE training";
 149   disp (tit);
 150   fflush (stdout);
 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);
 155   subplot (2, 2, 2);
 156   mesh (xi, yi, zm);
 157   title (tit);
 158   subplot(2, 2, 4);
 159   hold on
 160   contourf (xi, yi, zm, 20);
 161   plot (xs, ys, "+6");
 162   hold off
 163   demo_octgpr_pause (1, 4);
 164   close
 165
 166 endfunction
 167
 168 function demo_octgpr2 (ncnt = 50, npt = 500)
 169
 170   global prntfmt;
 171   npt = ncnt*ceil (npt/ncnt);
 172   U = rand (ncnt, 2);
 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);
 178
 179   [U, ur] = rbf_centers(X, ncnt);
 180
 181   fi = linspace (0, 2*pi, 20);
 182   ncolors = rows (colormap);
 183   hold on
 184   for i = 1:rows (U)
 185     xc = U(i,1) + ur(i) * cos (fi);
 186     yc = U(i,2) + ur(i) * sin (fi);
 187     line (xc, yc);
 188   endfor
 189   hold off
 190   demo_octgpr_pause (2, 2);
 191   close
 192
 193 endfunction
 194
 195 function demo_octgpr3 (ncnt = 100, nsamp = 1000)
 196
 197   global prntfmt;
 198   tit = "a peaked surface";
 199   disp (tit);
 200
 201   % create the mesh onto which to interpolate
 202   t = linspace (-3, 3, 50);
 203   [xi,yi] = meshgrid (t, t);
 204
 205   % evaluate
 206   zi = testfun1 (xi, yi);
 207   zimax = max (vec (zi)); zimin = min (vec (zi));
 208   subplot (2, 2, 1);
 209   mesh (xi, yi, zi);
 210   title (tit);
 211   subplot (2, 2, 3);
 212   contourf (xi, yi, zi, 20);
 213   demo_octgpr_pause (1, 1);
 214
 215   tit = sprintf ("sampled at %d random points, selected %d centers", nsamp, ncnt);
 216   disp (tit);
 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);
 222   xys = [xs ys];
 223
 224   % select centers using k-means
 225   xyc = rbf_centers (xys, ncnt);
 226   xc = xyc(:,1); yc = xyc(:,2);
 227
 228   subplot (2, 2, 2);
 229   plot3 (xs, ys, zs, ".+");
 230   title (tit);
 231   subplot (2, 2, 3);
 232   hold on
 233   plot (xs, ys, "+6");
 234   hold off
 235   subplot (2, 2, 4);
 236   hold on
 237   plot (xs, ys, "+");
 238   plot (xc, yc, "o2");
 239   hold off
 240   demo_octgpr_pause (1, 2);
 241
 242   tit = "PGP model with heuristic hypers";
 243   disp (tit);
 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);
 249   subplot (2, 2, 2);
 250   mesh (xi, yi, zm);
 251   title (tit);
 252   subplot(2, 2, 4);
 253   hold on
 254   contourf (xi, yi, zm, 20);
 255   plot (xs, ys, "+6");
 256   plot (xc, yc, "o5");
 257   hold off
 258   demo_octgpr_pause (1, 3);
 259
 260   tit = "PGP model with MLE training";
 261   disp (tit);
 262   fflush (stdout);
 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);
 267   subplot (2, 2, 2);
 268   mesh (xi, yi, zm);
 269   title (tit);
 270   subplot(2, 2, 4);
 271   hold on
 272   contourf (xi, yi, zm, 20);
 273   plot (xs, ys, "+6");
 274   plot (xc, yc, "o5");
 275   hold off
 276   demo_octgpr_pause (1, 4);
 277   close
 278
 279 endfunction