%%%%%
% S is an nx2 matrix of points in R2
% k is the number of clusters to find
% Sigma is the affinity variance to use
%
% Returns a vector of n cluster labels for each example
% with each label in [1:k]
%%%%
function [clusterLabels] = ngcluster( S, k, sigma )

    % Build affinity matrix, filling first with straight euclidean squared
    % distances

    n = size( S, 1 );
    for i=1:n
        for j=1:n
            A(i,j) = ( S(i,1) - S(j,1) )^2 + ( S(i,2) - S(j,2) )^2;
        end
        A( i, i ) = 0;
    end

    
    % now compute gaussian fall-off over
    % entire matrix, element by element
    A = exp( -A / (2*sigma^2) );

    % D's diagonal are the row sums from A
    D = diag( sum( A' ) );

    L = (D^(-1/2)) * A * (D^(-1/2));

    [eigVec, eigVal] = eig( L );

    
    X = eigVec(:,1:2);


    % Normalize each of X's rows
    Y = X ./ [ (sum( (X.^2)')').^(1/2) , (sum( (X.^2)')').^(1/2) ];

    % this is equivalent to:

    %for i=1:n 
    %    Y(i,1) = X(i,1) /  sqrt( X(i,1)^2 + X(i,2)^2 );
    %    Y(i,2) = X(i,2) /  sqrt( X(i,1)^2 + X(i,2)^2 );
    %end


    
    hold off
    subplot( 2, 1, 1 )
    plot( Y(:,1), Y(:,2), 'o' );

    % also plot unit circle, for reference
    hold on
    plot( cos( 0:0.1:2*pi), sin( 0:0.1:2*pi), 'r--' );


    % compute clusters in eigen RK using standard alg
    clusterLabels = clusterdata( Y, k );

    
    hold off
    subplot( 2, 1, 2 )

    colors = { 'ro', 'go', 'bo', 'po', 'ko' };

    % plot first few clusters
    plot( S( clusterLabels==1, 1), S( clusterLabels==1, 2), 'ro'  );
    hold on;
    if k>=2
        plot( S( clusterLabels==2, 1), S( clusterLabels==2, 2), 'go'  );
    end
    if k>=3
        plot( S( clusterLabels==3, 1), S( clusterLabels==3, 2), 'bo'  );
    end
    if k>=4
        plot( S( clusterLabels==4, 1), S( clusterLabels==4, 2), 'po'  );
    end

    hold off

