octave_packages/octgpr-1.2.0/rbf_centers.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} {[U, ur, iu]} = rbf_centers (@var{X}, @var{nu}, @var{theta})
  24 % Selects a given number of RBF centers based on Lloyd's clustering algorithm.
  25 %
  26 % @end deftypefn
  27 function [U, ur, iu] = rbf_centers (X, nu, theta)
  28
  29   if (nargin == 3)
  30     X *= diag (theta);
  31   elseif (nargin != 2)
  32     print_usage ();
  33   endif
  34
  35   % the D^2 weighting initialization
  36
  37   D = Inf;
  38   kk = 1:rows (X);
  39   cp = kk;
  40
  41   for i = 1:nu
  42     jj = sum (rand() * cp(end) < cp);
  43     k(i) = kk(jj);
  44     kk(jj) = [];
  45     U = X(k(i),:);
  46     D = min (D, pdist2_mw(X, U, 'ssq')');
  47     cp = cumsum (D(kk));
  48   endfor
  49
  50
  51   % now perform the k-means algorithm
  52
  53   U = X(k,:);
  54   D = pdist2_mw (U, X, 'ssq');
  55   [xx, j] = min (D);
  56
  57   it = 0;
  58   do
  59     for i = 1:columns (X)
  60       U(:,i) = accumarray (j.', X(:,i), [nu, 1]);
  61     endfor
  62     N = accumarray (j.', ones (1, length (j)), [nu, 1]);
  63     U = diag (N) \ U;
  64     i = find (all (U == 0, 2));
  65     U(i,:) = X(ceil (rand (1, length (i)) * rows (X)), :);
  66     j1 = j;
  67     D = pdist2_mw (U, X, 'ssq');
  68     [xx, j] = min (D);
  69     fprintf (stderr, "k-means iteration %d\r", ++it);
  70   until (all (j == j1))
  71   fprintf (stderr, "\n");
  72
  73   if (nargout > 2)
  74     iu = j;
  75   endif
  76
  77   if (nargout > 1)
  78     ur = zeros (nu, 1);
  79     for i = 1:nu
  80       ij = (j == i);
  81       ur(i) = sqrt (max (D(i,ij)));
  82     endfor
  83   endif
  84
  85   if (nargin == 3)
  86     U = dmult (U, 1./theta);
  87     if (any(theta == 0))
  88       U(:,theta == 0) = 0;
  89     endif
  90   endif
  91
  92 endfunction