%% Advanced Human Activity Clustering
% This script extends the example from ClusterAnalysis.m to evaluate other
% clustering methods. In particular, it addresses:
% * Hierarchical Clustering, Dendrograms & Cophenetic Correlation
% * Fuzzy C-means Clustering (Fuzzy Logic Toolbox)
% * Self-Organizing Maps (Neural Network Toolbox)
% * Evaluating & Comparing Clustering Solutions
% 
%  The dataset is courtesy:
%  Davide Anguita, Alessandro Ghio, Luca Oneto, Xavier Parra and Jorge L.
%  Reyes-Ortiz. Human Activity Recognition on Smartphones using a
%  Multiclass Hardware-Friendly Support Vector Machine. International
%  Workshop of Ambient Assisted Living (IWAAL 2012). Vitoria-Gasteiz,
%  Spain. Dec 2012

%% Load Data
% Load Human-Activity dataset for clustering.

load Data\HAData Xtrain features

% Perform PCA decomposition
[coeff,score,latent] = pca(Xtrain);

%% Initial Clustering Solution Baseline
% Compute initial K-means clustering solution for comparison

rng default
kidx = kmeans(Xtrain, 3, 'distance', 'correlation', 'replicates', 5);

%% Hierarchical Clustering, Cophenetic Corr. Coeff. and Dendrogram
% Hierarchical Clustering groups data over a variety of scales by creating
% a cluster tree or dendrogram. The tree is not a single set of clusters,
% but rather a multilevel hierarchy, where clusters at one level are joined
% as clusters at the next level. A variety of distance measures such as
% Euclidean, Mahalanobis, Jaccard etc. are available.
%
% Cophenetic correlation coefficient is typically used to evaluate which
% type of hierarchical clustering is more suitable for a particular type of
% data. It is a measure of how faithfully the tree represents the
% similarities/dissimilarities among observations.
% 
% Dendrogram is a graphical representation of the cluster tree created by
% hierarchical clustering.

Z = linkage(score(:,1:8), 'ward');

cpcc = cophenet(Z, pdist(score(:,1:8)));
disp('Cophenetic correlation coefficient: ')
disp(cpcc)

dendrogram(Z);
xlabel('Data point index')
ylabel ('Distance b/w data points')
title(['Cophenetic Correlation Coefficient: ' sprintf('%0.4f',cpcc)])

%%
% One can determine the appropriate number of clusters from the heights of
% the dendrogram. There seem to be two distinct clusters with perhaps 2-3
% sub-clusters each, so anywhere from 2 to 6 clusters may be appropriate

hidx = cluster(Z, 'maxclust', 3);
pairscatter(score(:,1:4), [], hidx)

%% 
% Evaluate clustering solution using silhouette values
clf
[sil,~] = silhouette(score(:,1:8), hidx);
title(sprintf('Average silhouette value: %0.2f', mean(sil)));

%% Fuzzy C-Means
% Fuzzy C-means (FCM) is a data clustering technique wherein each data
% point belongs to a cluster to some degree that is specified by a
% membership grade. The function _fcm_ returns a list of cluster centers
% and several membership grades for each data point.
% 
% _fcm_ iteratively minimizes the sum of distances of each data point to
% its cluster center weighted by that data point's membership grade.

[center,U] = fcm(score(:,1:8), 3);

clf
scatter(score(:,1), score(:,3), 3, U');
xlabel('PCA Score 1'); ylabel('PCA Score 3');
title('Fuzzy C-Means Cluster Membership Weight');
grid on;

% Generate hard clustering for evaluation
[~, fidx] = max(U);
fidx = fidx';

%% Neural Networks - Self Organizing Maps (SOMs)
% Neural Network Toolbox supports unsupervised learning with
% self-organizing maps (SOMs) and competitive layers.
% 
% SOMs learn to classify input vectors according to how they are grouped in
% the input space. SOMs learn both the distribution and topology of the
% input vectors they are trained on. One can make use of the interactive
% tools to setup, train and test a neural network. It is then possible to
% auto-generate the code for the purpose of automation. The code in this
% section has been auto-generated. 

% Create a Self-Organizing Map
dimension1 = 3;
dimension2 = 1;
net = selforgmap([dimension1 dimension2]);

% Train the Network
net.trainParam.showWindow = 1;
[net,tr] = train(net, Xtrain');

% Test the Network
nidx = net(Xtrain');
nidx = vec2ind(nidx)';

% Visualize results
pairscatter(score(:,1:4),[],nidx);

%% Gaussian Mixture Models (GMM)
% Gaussian mixture models are formed by combining multivariate normal
% density components. _gmdistribution.fit_ uses expectation maximization
% (EM) algorithm, which assigns posterior probabilities to each component
% density with respect to each observation. The posterior probabilities for
% each point indicate that each data point has some probability of
% belonging to each cluster.

gmobj = fitgmdist(score(:,1:8), 3, 'Replicates', 5);
gidx = cluster(gmobj, score(:,1:8));

clf
scatter(score(:,1), score(:,3), 3, posterior(gmobj, score(:,1:8)));
xlabel('PCA Score 1'); ylabel('PCA Score 3');
title('GMM Posterior Probabilities');
grid on;

figure
pairscatter(score(:,1:4), [], gidx);

%% Cluster Evaluation
% Compare clustering solutions using average silhouette value. This is one
% objective criterion. It may not guarantee the best clustering solution,
% which can be a subjective determination.

sols = [kidx hidx fidx nidx gidx];
methods = {'KMeans','Hierarchical','Fuzzy C-means','SOM','GMM'};
silAve = zeros(1, size(sols,2));
for i = 1:size(sols,2)
    subplot(2,3,i);
    %[silh,~] = silhouette(Xtrain, sols(:,i), 'correlation');
    [silh,~] = silhouette(score(:,1:8), sols(:,i));
    silAve(i) = mean(silh);
    title(sprintf('%s Ave Silhouette: %0.2f', methods{i}, silAve(i)));
end
subplot(2,3,6);
barh(silAve); xlabel('Average Silhouette Value');
set(gca,'YTick', 1:length(silAve), 'YTickLabel', methods);

%% Determining Number of Clusters
% The evalclusters function can be used to evaluate a clustering algorithm
% like kmeans for different number of clusters and calculate a metric that
% can be used to determine the optimal number of clusters. 
%
% Notice that here, for both kmeans and hierarchical clustering, the
% optimal number of clusters is 2. 

ev1 = evalclusters(score(:,1:8),'kmeans','CalinskiHarabasz','KList',2:6)
ev2 = evalclusters(score(:,1:6),'kmeans','DaviesBouldin','KList',2:6)

ev1 = evalclusters(score(:,1:8),'linkage','CalinskiHarabasz','KList',2:6)
ev2 = evalclusters(score(:,1:6),'linkage','DaviesBouldin','KList',2:6)


