function basis = kpca_calc(xs,kernel,d,kmataxis);%KPCA_CALC calculates a kernel PCA basis.%%   usage%      basis = kpca_calc(xs,kernel,d);%%   input%      xs       matrix of column vectors%      kernel   a chosen kernel, default = {'gaussian',1}%      d        number of eigenvectors (give for efficiency),%               default = size(xs,2)%      kmataxis is a figure handle where the kernel matrix will be%               plotted (default = 0 no plot)%%   output%      basis    struct containing the following entries%                   basis.V        eigenvectors%                   basis.Lambda   eigenvalues%                   basis.xs       used vectors%                   basis.kernel   used kernel%%   see also%      kpca_plot, kpca_map%%   STH * 12MAR2002if ~exist('kernel')|isempty(kernel), kernel = {'gaussian',1}; endif ~exist('d')|isempty(d), d = size(xs,2); endif ~exist('kmataxis')|isempty(kmataxis), kmataxis = 0; end% d can't be larger than the number of samplesif d>size(xs,2)  warning('d is larger than the number of samples, resetting d')  d = size(xs,2);endxsc = size(xs,2);   % column of xs% calculate the kernel matrixK = kpca_matrix(xs,xs,kernel);if kmataxis>0  cf = gcf;  figure(kmataxis)  imagesc(K)  figure(cf)end% center the kernel matrixsk = size(K,1);             % note, K is square matrixrowK = sum(K)/sk;           % the sums of the columnsallK = sum(K(:))/(sk*sk);   % the sum of all entriesK = K - repmat(rowK,[sk 1]) - repmat(rowK',[1 sk]) + repmat(allK,[sk sk]);% find the eigenvectors and eigenvaluesswitch 2 case 1  [V,Lambda] = jdqr(K/sk,d); case 2  opts.disp = 0;  [V,Lambda,flag] = eigs(K/sk,d,'LM',opts);  if flag    warning([mfilename ': not all eigenvalues converged'])  endend% we can not assume that the eigenvalues are sorted[dummy, ind] = sort(-diag(Lambda));Lambda = Lambda(ind,ind);V = V(:,ind);% due to numerical instabilities some eigenvalues might be negative% or smaller than eps, we want to ignore thosevalid = find(diag(Lambda)<2*eps);if length(valid)<1  % all eigenvalues are valid, keep d unchangedelse  % some are not valid  d = valid(1)-1;  warning([mfilename ': some eigenvalues of kernel matrix are less than eps'])endclear valid% cut off those eigenvalues and eigenvectorsV = V(:,1:d);Lambda = Lambda(1:d,1:d);% normalize the eigenvectors in feature spaceV = V*inv(sqrtm(sk*Lambda));% assign structbasis.V = V;basis.Lambda = Lambda;basis.xs = xs;basis.kernel = kernel;