MATLAB_DDC_Algorithm.m

function [ Clusters, Results ] = DDC_ver01_1_CAMS( varargin )
%DDC_VER01.1 Data Density Based Clustering
% Copyright R Hyde 2017
% Released under the GNU GPLver3.0
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/
% If you use this file please acknowledge the author and cite as a
% reference:
% Hyde, R.; Angelov, P., "Data density based clustering," Computational
% Intelligence (UKCI), 2014 14th UK Workshop on , vol., no., pp.1,7, 8-10 Sept. 2014
% doi: 10.1109/UKCI.2014.6930157
% Downloadable from: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6930157&isnumber=6930143
%
% ACCMIP revision: *Try clusters with no distance in any dimension set grid
% dimension instead of remaining as initial radii (Line 132)
%
% Data Density Based Clustering with Manual Radii
% Useage:
% [Clusters, Results]=DDC_ver01(DataIn, InitR, Merge, Verbose]
% Inputs:
% DataIn:   m x n array of data for clustering m rows of samples, data
%           should be normalised 0-1 or scaled appropriately
% InitR:    initial radii array with radius in each diension, if a single
%           value is provided the same value will be used for all.
% Verbose:  1 - plot output of progress, 0 - silent
% Merge:    Flag=1 to merge clusters if centre is within the ellipse of
%           another
% Outputs:
% Results:  array of data with cluster number appended as last column
% Clusters: array of cluster centre co-ords and radii

%% Start Main Function
if size(varargin,2)<4
    sprintf('Input error - not enough inputs.')
    return
end
DataIn=varargin{1}; % Read Array of Data
InitR=varargin{2}; % Read Initial Radius
Verbose=varargin{4}; % Read flag for plotting during the run
Merge=varargin{3}; % Read flag for merging final clusters
% defaults to use instead of initial radii, 1 grid in lat, lon, scaled O3 std dev in case of cluster radii = 0
% DefaultRadii = [getappdata(CatacombePaper,'GridDistance'),...
%     getappdata(CatacombePaper,  'RadO3')]; 

%% Initialise
if size(InitR,2)==1 % create equal radii across each domension if only one provided
    if Verbose==1
    sprintf('Using equal radii')
    end
    InitR=ones(size(DataIn,2),1)*InitR;
elseif size(InitR,2)<size(DataIn,2)
    sprintf('Radii not single or equal to data dimensions')
    return
end
NumClusters=0; % initial number of clusters
Results=zeros(size(DataIn,1),size(DataIn,2)+1); % initiate empty array

%% ### DDC routine ###
Glob_Mean=mean(DataIn,1); % array of means of data dim
Glob_Scalar=sum(sum((DataIn.*DataIn),2),1)/size(DataIn,1); % array of scalar products for each data dim

while size(DataIn,1)>0
    % size(DataIn,1) % uncomment to trace remaining data
    NumClusters=NumClusters+1;
    Clusters.Rad(NumClusters,:)=InitR;
    %% Find Cluster Centre
    Glob_Mean=mean(DataIn,1); % array of means of data dim
    Glob_Scalar=sum(sum((DataIn.*DataIn),2),1)/size(DataIn,1); % array of scalar products for each data dim
  % full calculations
%     GDensity=1./(1+(pdist2(DataIn,Glob_Mean,'euclidean').^2)+Glob_Scalar-(sum(Glob_Mean.^2))); % calculate global densities
%     [~, CentreIndex]=max(GDensity); % find index of max densest point
    % slim calculations
    GDensity=pdist2(DataIn,Glob_Mean,'euclidean').^2 + Glob_Scalar - sum(Glob_Mean.^2); % calculate global densities
    [~, CentreIndex]=min(GDensity); % find index of max densest point
  
    %% Find points belonging to cluster
    Include=bsxfun(@minus,DataIn,DataIn(CentreIndex,:)).^2; % sum square of distances from centre
    RadSq=Clusters.Rad(NumClusters,:).^2; % square radii
    Include=sum(bsxfun(@rdivide,Include,RadSq),2); % divide by radii and add terms
    Include=find(Include<1);
    
    %% Remove outliers >3sigma
    Dist=pdist2(DataIn(Include,:),DataIn(CentreIndex,:)); % distances to all potential members
    Include=Include(abs(Dist - mean(Dist) <= 3*std(Dist))==1,:); % keep only indices of samples with 3 sigma

    %% Move cluster centre to local densest point
    LocMean=mean(DataIn(Include,:),1);
    LocScalar=sum((DataIn(Include,:).^2),2)/size(Include,1); % array of scalar products of data dims
    % full calculations
%     LocDens=1./(1+(pdist2(DataIn(Include,:),LocMean,'euclidean').^2)+LocScalar-(sum(LocMean.^2))); % calculate local densities
%     [~,CentreIndex]=max(LocDens);
    % slim calculations
    LocDens=pdist2(DataIn(Include,:),LocMean,'euclidean').^2 + LocScalar - sum(LocMean.^2); % calculate local densities
    [~,CentreIndex]=min(LocDens);
    CentreIndex=Include(CentreIndex);
    Clusters.Centre(NumClusters,:)=DataIn(CentreIndex,:); % assign cluster centre
    
    %% Assign data to new centre
    Include=bsxfun(@minus,DataIn,Clusters.Centre(NumClusters,:)).^2; % sum square of distances from centre
    RadSq=Clusters.Rad(NumClusters,:).^2; % square radii
    Include=sum(bsxfun(@rdivide,Include,RadSq),2); % divide by radii and add terms
    Include=find(Include<1);

    %% Remove outliers >3sigma
    Dist=pdist2(Clusters.Centre(NumClusters,:),DataIn(Include,:)); % distances to all potential members
    Include=Include(abs(Dist - mean(Dist) <= 3*std(Dist))==1,:); % keep only indices of samples with 3 sigma

    %% Update radii to maximum distances
    for idx=1:size(DataIn,2)
        value01=pdist2(DataIn(Include,idx),Clusters.Centre(NumClusters,idx),'Euclidean');
        if max(value01)>0
            Clusters.Rad(NumClusters,idx)=max(value01);
        end
    end
    
    %% Assign data to cluster based on new radii
    Include=bsxfun(@minus,DataIn,Clusters.Centre(NumClusters,:)).^2; % sum square of distances from centre
    RadSq=Clusters.Rad(NumClusters,:).^2; % square radii
    Include=sum(bsxfun(@rdivide,Include,RadSq),2); % divide by radii and add terms
    Include=find(Include<1);
 
    %% Remove outliers >3sigma
    Dist=pdist2(Clusters.Centre(NumClusters,:),DataIn(Include,:)); % distances to all potential members
    Include=Include(abs(Dist - mean(Dist) <= 3*std(Dist))==1,:); % keep only indices of samples with 3 sigma
    
    %% Update radii to maximum distances
    
    for idx=1:size(DataIn,2)
        value01=pdist2(DataIn(Include,idx),Clusters.Centre(NumClusters,idx),'Euclidean');
        if max(value01)>0
            Clusters.Rad(NumClusters,idx)=max(value01);
        else
%             Clusters.Rad(NumClusters,idx)=DefaultRadii(idx);
        end
    end
 
    %% Plot
    if Verbose==1
        hold off;scatter(DataIn(:,1),DataIn(:,2));hold on
        scatter(DataIn(CentreIndex,1),DataIn(CentreIndex,2),'r')
        scatter(DataIn(Include,1),DataIn(Include,2),'g');
        scatter(Clusters.Centre(NumClusters,1),Clusters.Centre(NumClusters,2),'*','r')
        title(sprintf('Clustered: %i, Remaining: %i',size(Results,1)-size(DataIn,1), size(DataIn,1)))
        axis([0 1 0 1])
        drawnow
        for zz=1:size(Clusters.Centre,1)
            rectangle('Position',[Clusters.Centre(zz,1)-Clusters.Rad(zz,1), Clusters.Centre(zz,2)-Clusters.Rad(zz,2), 2*Clusters.Rad(zz,1), 2*Clusters.Rad(zz,2)],'Curvature',[1,1])
        end
    end
    %% Assign data to final clusters
    StartIdx=find(all(Results==0,2),1,'first');
    EndIdx=StartIdx+size(Include,1)-1;
    Results(StartIdx:EndIdx,:)=[DataIn(Include,:),ones(size(Include,1),1)*NumClusters];
    DataIn(Include,:)=[]; % remove clustered data
end

%% Merge clusters if centre is within another cluster
if Merge==1;
MergeAny=1;
    while MergeAny==1
        if Verbose==1
            figure(2)
            clf
            for zz=1:size(Clusters.Centre,1)
                rectangle('Position',[Clusters.Centre(zz,1)-Clusters.Rad(zz,1),...
                    Clusters.Centre(zz,2)-Clusters.Rad(zz,2), 2*Clusters.Rad(zz,1),...
                    2*Clusters.Rad(zz,2)],'Curvature',[1,1])
            end
            hold on
            scatter(Clusters.Centre(:,1),Clusters.Centre(:,2),'*','r')
            drawnow
        end

        MergeAny=0;
        Merges=[];
        % for each cluster & find if cluster centre is within other clusters
        for idx1=1:size(Clusters.Centre,1); 
            InEll=bsxfun(@minus,Clusters.Centre,Clusters.Centre(idx1,1:end)).^2;
            InEll=sum(bsxfun(@rdivide,InEll,Clusters.Rad(idx1,:).^2),2); % divide by rad^2 & add
            InEll=(InEll<1);
            Merges(idx1,:)=InEll.';
        end
        Merges(logical(eye(size(Merges))))=0;
        % Merge clusters
        for idx=1:size(Clusters.Centre,1)
            [~,idx1]=find(Merges(idx,:),1);
            Results(ismember(Results(:,end),idx1),end)=idx;
            if idx1
                MergeAny=1;
            end
        end
        %% renumber clusters
        [C,~,ic]=unique(Results(:,end));
        C=1:size(C,1);
        Results(:,end)=C(ic);
        %% Re-create cluster data
        Clusters.Centre=[];
        Clusters.Rad=[];
        for idx1=1:max(Results(:,end))
            Clusters.Centre(idx1,:)=mean(Results(Results(:,3)==idx1,1:end-1),1);
            for idx2=1:size(Results,2)-1
                value01=pdist2(Results(Results(:,3)==idx1,idx2),Clusters.Centre(idx1,idx2),'Euclidean');
                if max(value01)>0
                    Clusters.Rad(idx1,idx2)=max(value01);
                else
                    Clusters.Rad(idx1,idx2)=0;
                end
            end
        end

    end
end

end % end function