function [patterns, targets, label] = spectral_k_means(train_patterns, train_targets, params, plot_on)

%Reduce the number of data points using the spectral k-means algorithm
%The k largest eigenvectors of the kernel matrix of examples are clustered
%using k-means
%
%Inputs:
%	train_patterns	- Input patterns
%	train_targets	- Input targets
%	params			- [Number of output data points, Kernel type, kernel parameter, cluster_type]
%                     Kernel can be one of: Gauss, RBF (Same as Gauss), Poly, Sigmoid, or Linear
%                     The kernel parameters are:
%                       RBF kernel  - Gaussian width (One parameter)
%                       Poly kernel - Polynomial degree
%                       Sigmoid     - The slope and constant of the sigmoid (in the format [1 2], with no separating commas)
%					    Linear		- None needed                   
%                     Clustering type can be Multicut (Meila-Shi) or NJW (Ng, Jordan, Weiss)
%   plot_on         - Unused
%
%Outputs
%	patterns		- New patterns
%	targets			- New targets
%	label			- The labels given for each of the original patterns

if (nargin < 4),
    plot_on = 0;
end

%Get parameters
[Nmu, kernel, kernel_params, clustering_type] = process_params(params);

max_iter   = 1000;
[Ndim, Nf] = size(train_patterns);
label      = zeros(1,Nf);

%Kernelize the training patterns
y       = kernelize(kernel, kernel_params, train_patterns, train_patterns);
D       = diag(sum(y));

switch lower(clustering_type)
    case 'multicut'
        y       = inv(D)*y;
        [v, d]  = eig(y);
        kernel_patterns = v(end-Nmu+1:end, :);
    case 'njw'
        sqD             = sqrtm(D);
        L               = sqD*y*sqD;
        [v, d]          = eig(L);
        kernel_patterns = v(end-Nmu+1:end, :);
        kernel_patterns = kernel_patterns ./ (sum(kernel_patterns'.^2)'*ones(1, Nf));
    otherwise
        error('Unknown clustering type');
end

[p, t, label] = k_means(kernel_patterns, train_targets, Nmu, []);

Ul      = unique(label);
patterns= zeros(Ndim, length(Ul));
targets = zeros(1,    length(Ul));
for i = 1:length(Ul)
    patterns(:, i) = mean(train_patterns(:, find(label == Ul(i)))')';
    targets(:, i)  = mean(train_targets(find(label == Ul(i)))')';
end    





%END kernel k-means

function y = kernelize(kernel, kernel_params, patterns, base_patterns)

Nf  = size(patterns, 2);
Nb  = size(base_patterns, 2);
y	= zeros(Nb, Nf);

switch kernel,
case {'Gauss','RBF'},
   for i = 1:Nf,
      y(:,i)    = exp(-sum((base_patterns-patterns(:,i)*ones(1,Nb)).^2)'/(2*kernel_params^2));
   end
case {'Poly', 'Linear'}
   if strcmp(kernel, 'Linear')
      kernel_params = 1;
   end
   
   for i = 1:Nf,
      y(:,i) = (base_patterns'*patterns(:,i) + 1).^kernel_params;
   end
case 'Sigmoid'
    kernel_params = str2num(kernel_params);
    
    if (length(kernel_params) ~= 2)
        error('This kernel needs two parameters to operate!')
    end
    
   for i = 1:Nf,
      y(:,i) = tanh(base_patterns'*patterns(:,i)*kernel_params(1)+kernel_params(2));
   end
otherwise
   error('Unknown kernel. Can be Gauss, Linear, Poly, or Sigmoid.')
end