%% K-means Algorithm for Angular Data
%
% Apply standard clustering algorithm using custom distance metric for
% angular (circular) data.
%
% Default distance metric: 1 - cos( theta_i - theta_j )
% Default similarity metric: cos( theta_i - theta_j )
%
% Standard K-means:
%
% 1) random initialization of cluster "means" to randomly-selected data
% points
%
% 2)compute data point associations based on "nearest" cluster "mean"
%
% 3)compute mean of candidate clusters
%
% 4)repeat (2) and (3) until convergence
%
% References:
%1]Clustangles: An Open Library for Clustering Angular Data
%  Karen Sargsyan*?, Yun Hao Hua?, and Carmay Lim*??
%  Institute of Biomedical Sciences, Academia Sinica, Taipei 115, Taiwan
%  Department of Chemistry, National Tsing Hua University, Hsinchu 300, Taiwan


%% Input Arguments

% K - Number of clusters
%     type: positive integer

% data - Angle (radians) data points. Format [sin(theta_1), cos(theta_1); ...
%        sin(theta_N), cos(theta_N)];
%        type: Nx2 matrix

% isLines - Optional special toggle for "line-orientation" clustering,
%           since each line data point corresponds to 2 possible
%           orientations requiring special handling of averaging and
%           distance measure.
%           type: boolean

% dist - Optional custom distance metric. Format: inputs=>{data(x,y), mean},
%        outputs=>{distance measure}.
%        type: function_handle

% tol - Optional convergence tolerance of cluster means. K-means iterations
%       terminate if change in "cluster-means" < angular-tolerance.
%       type: positive scalar

% maxIter - Optional maximum number of Kmeans iterations.
%           type: positive integer


%% Output Arguments

% assoc - Data point cluster associations/assignments. 
%         type: Nx1 vector


%% K-means Algorithm for Angular Data

function assoc = angular_kmeans(K, data, varargin)

debugEnable = false;

% Check input parameters
[K, data, nPts, tol, maxIter, dist, isLines] = input_parser(K, data, varargin{:});

% K-means Algorithm
prevMeans = nan(K, 2);

% Random initialization of means
currMeans = data(randperm(nPts,K), :);

niter = 1;
while( ~all(abs(acos(diag(currMeans*prevMeans')))<tol) && niter <= maxIter )

    prevMeans = currMeans;
    
    % Compute cluster associations for all data-points
    if isLines,
        simMat = abs(data*currMeans');
        [~, assoc] = max( simMat, [], 2);
    elseif isempty(dist),
        simMat = data*currMeans';
        [~, assoc] = max( simMat, [], 2);
    else
        distMat = nan(nPts, K);
        for i = 1:nPts,
            for j = 1:K,
                distMat = dist(data(i,:), currMeans(j,:));
            end
        end
        [~, assoc] = min( distMat, [], 2);
    end    
    
    % For "line-orientation"-only, flip the line vectors which achieved a
    % higher similarity in the flipped orientation prior to computing
    % Euclidean means
    if isLines,
       
        for i = 1:nPts,
            if data(i,:)*currMeans(assoc(i),:)' < 0 
                data(i,:) = -data(i,:);
            end
        end
        
    end
    
    % Compute cluster "means"               
    for i = 1:K,
        ave = mean(data(assoc==i, :));
        ang = atan2( ave(1), ave(2) );
        currMeans(i,:) = [sin(ang), cos(ang)];
    end 
    
    niter = niter + 1;
    
end

% Plot data points and different highlighting for each cluster
if debugEnable
    
   figure;
   ax1 = gca;
      
   colors = jet(K);
   for c = 1:K,
              
      hold on; scatter(ax1, data(assoc==c, 2), data(assoc==c, 1), 8, colors(c,:));
      plot(currMeans(c,2), currMeans(c,1), '*', 'color', colors(c,:));
      
   end    
    
end


end


%% Input Parser

function [K, data, nPts, tol, maxIter, dist, isLines] = input_parser(K, data, varargin)

% parse input arguments
p = inputParser;
p.addRequired('K', @(x) ...
              validateattributes(x,{'numeric'},{'real','integer', 'positive'}));
p.addRequired('data', @(x) ...
              validateattributes(x,{'numeric'},{'real', '2d', 'nonnan', 'finite', 'ncols', 2}));
p.addParameter('tol', 1e-3, @(x) ...
              validateattributes(x,{'numeric'},{'real', 'scalar', 'positive', 'finite', 'nonnan'}));          
p.addParameter('maxIter', 100, @(x) ...
              validateattributes(x,{'numeric'},{'real','integer', 'positive'}));          
p.addParameter('dist', [], @(x) isa(x, 'function_handle'));  
p.addParameter('isLines', false, @(x) ...
              validateattributes(x,{'logical'},{'real','nonnegative'}));    
p.parse(K, data, varargin{:});
inputs       = p.Results;
K            = inputs.K;
data         = inputs.data;
tol          = inputs.tol;
maxIter      = inputs.maxIter;
dist         = inputs.dist;
isLines      = inputs.isLines;

nPts = size(data, 1);

end %end of input parser












