X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Foctgpr-1.2.0%2Fdemo_octgpr.m;fp=octave_packages%2Foctgpr-1.2.0%2Fdemo_octgpr.m;h=3482230f27a273cf47dd1ed45bc3c4b578055459;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git

diff --git a/octave_packages/octgpr-1.2.0/demo_octgpr.m b/octave_packages/octgpr-1.2.0/demo_octgpr.m
new file mode 100644
index 0000000..3482230
--- /dev/null
+++ b/octave_packages/octgpr-1.2.0/demo_octgpr.m
@@ -0,0 +1,279 @@
+% Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
+% 
+% Author: Jaroslav Hajek <highegg@gmail.com>
+% 
+% This file is part of OctGPR.
+% 
+% OctGPR 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.
+% 
+% This program 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 this software; see the file COPYING.  If not, see
+% <http://www.gnu.org/licenses/>.
+% 
+
+% -*- texinfo -*- 
+% @deftypefn {Function File} demo_octgpr (1, nsamp = 150)
+% @deftypefnx {Function File} demo_octgpr (2, ncnt = 20, npt = 500)
+% @deftypefnx {Function File} demo_octgpr (3, ncnt = 50, nsamp = 500)
+% OctGPR package demo function.
+% First argument selects available demos:
+%
+% @itemize
+% @item 1. GPR regression demo @*
+% A function is sampled (with small noise), then reconstructed using GPR
+% regression.  @var{nsamp} specifies the number of samples.
+% @seealso{gpr_train, gpr_predict}
+% @item 2. RBF centers selection demo @*
+% Radial basis centers are selected amongst random points.
+% @var{ncnt} specifies number of centers, @var{npt} number of points.
+% @item 2. PGP regression demo @*
+% A function is densely sampled (with small noise), 
+% radial basis centers are selected, then the function is reconstructed 
+% using PGP regression.  @var{nsamp} specifies the number of samples,
+% @var{ncnt} specifies number of centers.
+% @end itemize
+% @end deftypefn
+function demo_octgpr (number, varargin)
+  
+  global prntfmt = '';
+
+  if (nargin < 1)
+    print_usage ();
+  elseif (ischar (number))
+    prntfmt = number;
+  elseif (isscalar (number))
+    figure ();
+    if (! isempty (prntfmt))
+      figure (gcf, "visible", "off");
+    endif
+
+    switch (number)
+    case 1
+      demo_octgpr1 (varargin{:})
+    case 2
+      demo_octgpr2 (varargin{:})
+    case 3
+      demo_octgpr3 (varargin{:})
+    otherwise
+      error ("demo_octgpr: invalid demo number")
+    endswitch
+  else
+    print_usage ();
+  endif
+
+endfunction
+
+function demo_octgpr_pause (idemo, iplot)
+  global prntfmt;
+  if (isempty (prntfmt))
+    pause;
+  else
+    print (sprintf (prntfmt, idemo, iplot));
+    fflush (stdout);
+  endif
+endfunction
+
+% define the test function (the well-known matlab "peaks" plus some sines)
+function z = testfun1 (x, y)
+  z = 4 + 3 * (1-x).^2 .* exp(-(x.^2) - (y+1).^2) ...
+      + 10 * (x/5 - x.^3 - y.^5) .* exp(-x.^2 - y.^2) ...
+      - 1/3 * exp(-(x+1).^2 - y.^2) ...
+      + 2*sin (x + y + 1e-1*x.*y);
+endfunction
+
+function demo_octgpr1 (nsamp = 150)
+  global prntfmt;
+  tit = "a peaked surface";
+  disp (tit);
+
+  % create the mesh onto which to interpolate
+  t = linspace (-3, 3, 50);
+  [xi,yi] = meshgrid (t, t);
+
+  % evaluate
+  zi = testfun1 (xi, yi);
+  zimax = max (vec (zi)); zimin = min (vec (zi));
+  subplot (2, 2, 1);
+  mesh (xi, yi, zi);
+  title (tit);
+  subplot (2, 2, 3);
+  contourf (xi, yi, zi, 20);
+  demo_octgpr_pause (1, 1);
+
+  tit = sprintf ("sampled at %d random points", nsamp);
+  disp (tit);
+  % create random samples
+  xs = rand (nsamp,1); ys = rand (nsamp,1);
+  xs = 6*xs-3; ys = 6*ys - 3;
+  % evaluate at random samples
+  zs = testfun1 (xs, ys);
+  xys = [xs ys];
+
+  subplot (2, 2, 2);
+  plot3 (xs, ys, zs, ".+");
+  title (tit);
+  subplot (2, 2, 3);
+  hold on
+  plot (xs, ys, "+6");
+  hold off
+  subplot (2, 2, 4);
+  plot (xs, ys, ".+");
+  demo_octgpr_pause (1, 2);
+
+  tit = "GPR model with heuristic hypers";
+  disp (tit);
+  ths = 1 ./ std (xys);
+  GPM = gpr_train (xys, zs, ths, 1e-5);
+  zm = gpr_predict (GPM, [vec(xi) vec(yi)]);
+  zm = reshape (zm, size(zi));
+  zm = min (zm, zimax); zm = max (zm, zimin);
+  subplot (2, 2, 2);
+  mesh (xi, yi, zm);
+  title (tit);
+  subplot(2, 2, 4);
+  hold on
+  contourf (xi, yi, zm, 20);
+  plot (xs, ys, "+6");
+  hold off
+  demo_octgpr_pause (1, 3);
+
+  tit = "GPR model with MLE training";
+  disp (tit);
+  fflush (stdout);
+  GPM = gpr_train (xys, zs, ths, 1e-5, {"tol", 1e-5, "maxev", 400, "numin", 1e-8});
+  zm = gpr_predict (GPM, [vec(xi) vec(yi)]);
+  zm = reshape (zm, size (zi));
+  zm = min (zm, zimax); zm = max (zm, zimin);
+  subplot (2, 2, 2);
+  mesh (xi, yi, zm);
+  title (tit);
+  subplot(2, 2, 4);
+  hold on
+  contourf (xi, yi, zm, 20);
+  plot (xs, ys, "+6");
+  hold off
+  demo_octgpr_pause (1, 4);
+  close
+
+endfunction
+
+function demo_octgpr2 (ncnt = 50, npt = 500)
+
+  global prntfmt;
+  npt = ncnt*ceil (npt/ncnt);
+  U = rand (ncnt, 2);
+  cs = min (pdist2_mw (U, 2) + diag (Inf (ncnt, 1)));
+  X = repmat (U, npt/ncnt, 1) + repmat (cs', npt/ncnt, 2) .* randn (npt, 2);
+  disp ("slightly clustered random points")
+  plot (X(:,1), X(:,2), "+");
+  demo_octgpr_pause (2, 1);
+
+  [U, ur] = rbf_centers(X, ncnt);
+
+  fi = linspace (0, 2*pi, 20);
+  ncolors = rows (colormap);
+  hold on
+  for i = 1:rows (U)
+    xc = U(i,1) + ur(i) * cos (fi);
+    yc = U(i,2) + ur(i) * sin (fi);
+    line (xc, yc);
+  endfor
+  hold off
+  demo_octgpr_pause (2, 2);
+  close
+
+endfunction
+
+function demo_octgpr3 (ncnt = 100, nsamp = 1000)
+
+  global prntfmt;
+  tit = "a peaked surface";
+  disp (tit);
+
+  % create the mesh onto which to interpolate
+  t = linspace (-3, 3, 50);
+  [xi,yi] = meshgrid (t, t);
+
+  % evaluate
+  zi = testfun1 (xi, yi);
+  zimax = max (vec (zi)); zimin = min (vec (zi));
+  subplot (2, 2, 1);
+  mesh (xi, yi, zi);
+  title (tit);
+  subplot (2, 2, 3);
+  contourf (xi, yi, zi, 20);
+  demo_octgpr_pause (1, 1);
+
+  tit = sprintf ("sampled at %d random points, selected %d centers", nsamp, ncnt);
+  disp (tit);
+  % create random samples
+  xs = rand (nsamp,1); ys = rand (nsamp,1);
+  xs = 6*xs-3; ys = 6*ys - 3;
+  % evaluate at random samples
+  zs = testfun1 (xs, ys);
+  xys = [xs ys];
+
+  % select centers using k-means
+  xyc = rbf_centers (xys, ncnt);
+  xc = xyc(:,1); yc = xyc(:,2);
+
+  subplot (2, 2, 2);
+  plot3 (xs, ys, zs, ".+");
+  title (tit);
+  subplot (2, 2, 3);
+  hold on
+  plot (xs, ys, "+6");
+  hold off
+  subplot (2, 2, 4);
+  hold on
+  plot (xs, ys, "+");
+  plot (xc, yc, "o2");
+  hold off
+  demo_octgpr_pause (1, 2);
+
+  tit = "PGP model with heuristic hypers";
+  disp (tit);
+  ths = 1 ./ std (xyc);
+  GPM = pgp_train (xys, xyc, zs, ths, 1e-5);
+  zm = pgp_predict (GPM, [vec(xi) vec(yi)]);
+  zm = reshape (zm, size(zi));
+  zm = min (zm, zimax); zm = max (zm, zimin);
+  subplot (2, 2, 2);
+  mesh (xi, yi, zm);
+  title (tit);
+  subplot(2, 2, 4);
+  hold on
+  contourf (xi, yi, zm, 20);
+  plot (xs, ys, "+6");
+  plot (xc, yc, "o5");
+  hold off
+  demo_octgpr_pause (1, 3);
+
+  tit = "PGP model with MLE training";
+  disp (tit);
+  fflush (stdout);
+  GPM = pgp_train (xys, xyc, zs, ths, 1e-3, {"tol", 1e-5, "maxev", 400});
+  zm = pgp_predict (GPM, [vec(xi) vec(yi)]);
+  zm = reshape (zm, size (zi));
+  zm = min (zm, zimax); zm = max (zm, zimin);
+  subplot (2, 2, 2);
+  mesh (xi, yi, zm);
+  title (tit);
+  subplot(2, 2, 4);
+  hold on
+  contourf (xi, yi, zm, 20);
+  plot (xs, ys, "+6");
+  plot (xc, yc, "o5");
+  hold off
+  demo_octgpr_pause (1, 4);
+  close
+
+endfunction