%Data = n x m matrix in which each row (1:n) contains observations for a single
%neuron for m total independent variables

%Step 1: calculate distances between data points in m-dimensional space
%using Euclidean distance measurement (other measurement options available,
%see Matlab Help search: pdist)
Data_dist = pdist(Data,'euclidean');

%Step 2: determine clustering based on distances using Ward's method
Data_link = linkage(Data_dist,'ward');

%Step 3: create a hierarchical cluster tree where N is the number of
%neurons (length of rows in n x m Data matrix)
Data_dend = dendrogram(Data_link,N);

%Step 4: generate vectors containing the row index for neurons in each
%cluster, e.g. cluster1 = [14 3 7 10 25...] based on dendrogram

%Step 5: make a vector with cluster asignments for each neuron:
groups = zeros(1,N);
groups(cluster1) = 1;
groups(cluster2) = 2;
groups(cluster3) = 3;
%etc depending on number of clusters identified in dendrogram

%Step 6: verify optimal number of clusters using combination of PCA and GMM
[~,score] = pca(Data); %output is eigenvalue "scores" for each neuron

GMModels = cell(3,1); %Preallocation assuming testing of 3 GMMs
options = statset('MaxIter',1000);
rng(1); %For reproducibility

for j = 1:3 %Displays mean cluster values for 3 GMMs with 1-3 putative clusters
    GMModels{j} = fitgmdist(score,j,'Options',options);
    fprintf('\n GM Mean for %i Component(s)\n',j)
    Mu = GMModels{j}.mu %displays mean values for each GMM
end

figure
    subplot(1,3,1), gscatter(score(:,1),score(:,2),groups); %plots clusters based first 3 PCs
    subplot(1,3,2), gscatter(score(:,2),score(:,3),groups);
    subplot(1,3,3), gscatter(score(:,1),score(:,3),groups);
  
%Step 7: determine minimum values indicating optimal number of clusters
[value,optimalclusternum] = min([GMModels{1}.AIC GMModels{2}.AIC GMModels{3}.AIC])
[value,optimalclusternum] = min([GMModels{1}.BIC GMModels{2}.BIC GMModels{3}.BIC])
[value,optimalclusternum] = min([GMModels{1}.NegativeLogLikelihood GMModels{2}.NegativeLogLikelihood GMModels{3}.NegativeLogLikelihood])

%Step 8: alternative - evaluate clusters using GMM approach with Calinski/Harabasz
%criterion and testing 1-6 putative clusters in Data        
rng('default');  %For reproducibility
clustereval = evalclusters(Data,'gmdistribution','CalinskiHarabasz','KList',[1:6]) 