%% clust.m
%%
%% CSS 581, Introduction to Machine Learning, Winter 2014
%% Computing & Software Systems, University of Washington Bothell
%% written by J. Jeffry Howbert, 2012-02-18
%%
%% k-means and hierarchical clustering on mixture of synthetic Gaussians

clear all;
close all;

load synthGaussMix

tic;

kMin = 2;
kMax = 8;

%% ***** for Exercises 6, you only need to change parameters in this block
distanceType = 'sqEuclidean';   % for kmeans: allowed values: 'sqEuclidean', ...
                                %   'cityblock', 'cosine', 'correlation'
linkType = 'average';            % for linkage: allowed values: 'single', ...
                                %   'complete', 'average'
doKmeans = true;                % if false, do hierarchical clustering instead
%% ***** end parameter block

nSamp = size( data, 1 );

for k = kMin : kMax
    fprintf( 1, '**** k = %d *****\n', k ); 
    % for kmeans, initial centroids chosen randomly, so important to set RandStream
    randSeed = 1;
    RandStream.setDefaultStream( RandStream( 'mt19937ar', 'Seed', randSeed ) );

    if ( doKmeans )
        % do k-means clustering
        clusIdx = kmeans( data, k, 'distance', distanceType, ...
            'start', 'sample' );
    else
        % do hierarchical clustering
        Z = linkage( data, linkType );
        % convert hierarchical clustering to partitional clustering for analysis
        clusIdx = cluster( Z, 'maxclust', k );
    end

    %% ***** unsupervised measures of cluster quality *****
    % SSE: sum of squared errors (error = Euclidean distance from
    %   assigned centroid).  Small value indicates higher quality.
    SSE = 0.0;
    for clus = 1 : k
        dataClus = data( clusIdx == clus, : );
        centroid = mean( dataClus, 1 );
        error = pdist2( dataClus, centroid, 'euclidean' );
        SSE = SSE + sum( error .^ 2 );
    end
    fprintf( 1, '   sum of squared error   = %8.4f\n', SSE );
    % Correlation of proximity and incidence matrices.  High absolute value
    %   (in range 0 to 1) indicates higher quality.
    dissim = squareform( pdist( data ) );   % use Euclidean distance as measure of dissimilarity
    incidence = zeros( nSamp, nSamp );
    for samp1 = 1 : nSamp
        for samp2 = 1 : nSamp
            if ( clusIdx( samp1 ) == clusIdx( samp2 ) )
                incidence( samp1, samp2 ) = 1.0;
            end
        end
    end
    R = corrcoef( dissim( : ), incidence( : ) );
    absCorr = abs( R( 1, 2 ) );
    fprintf( 1, '   correlation            = %6.4f\n', absCorr );

    % Print distribution of datapoints across clusters
    fprintf( 1, '   datapoints per cluster =' );
    for clus = 1 : k
        fprintf( 1, '%5d', sum( clusIdx == clus ) );
    end
    fprintf( 1, '\n' );    
    
end

toc;
