 classdef (Abstract) ClusterCriterion < classreg.learning.internal.DisallowVectorOps...
                           & matlab.mixin.CustomDisplay
% ClusterCriterion Clustering evaluation object.
%   ClusterCriterion is an abstract class used for clustering evaluation.
%   You cannot create an instance of this class directly. Create an
%   instance of a derived class such as SilhouetteEvaluation,
%   GapEvaluation, CalinskiHarabaszEvaluation or DaviesBouldinEvaluation by
%   calling the function EVALCLUSTERS
%
%   ClusterCriterion properties:
%      X                     - Data used to create this object.
%      NumObservations       - Number of observations in X.
%      OptimalK              - The optimal number of clusters suggested.
%      OptimalY              - The optimal clustering solution.
%      InspectedK            - List of the number of clusters inspected.
%      CriterionName         - Name of the criterion.
%      CriterionValues       - The criterion values for each number of clusters.
%      ClusteringFunction    - Clustering algorithm used.
%      Missing               - A vector indicating which observations are missing.
%
%   ClusterCriterion methods
%      compact               - Compact this object.
%      plot                  - Plot the criterion values.
%      addK                  - Add a new list of number of clusters for
%                              evaluation.

%   Copyright 2013-2023 The MathWorks, Inc.


properties(GetAccess=public, SetAccess=protected) 
%OptimalY The suggested clustering solution.
%   The OptimalY property is a column vector saving the clustering solution
%   corresponding to the value of OptimalK. It will be empty if the
%   clustering solutions are provided in the input.
    OptimalY = [];

%X matrix of data used to create this evaluation object.
%   The X property is a numeric matrix of data with one row per observation
%   and one column per variable.
    X;
    
%InspectedK The list of number of clusters that have been inspected.
    %   The InspectedK property is a vector of integers indicating the
    %   number of clusters inspected.
    InspectedK;
    
%CriterionValues The criterion values.
%   The CriterionValues property saves the criterion values corresponding
%   to each number of clusters in the InspectedK property.
    CriterionValues;
    
%CriterionName Name of the criterion. 
%   The CriterionName property is a string indicating the name of the
%   criterion requested for evaluation.
    CriterionName;
    
%OptimalK The suggested number of clusters.
%   The OptimalK property is an integer indicating the suggested number
%   of clusters after comparing the values of the requested criterion.
%   The value of OptimalK is estimated based on the following:
%   'CalinskiHarabasz'    The number of clusters giving the
%                         maximum CalinskiHarabasz value.
%   'silhouette'          The number of clusters giving the maximum
%                         silhouette value.
%   'gap'                 The number of clusters giving the first local
%                         maximum gap value or the global maximum gap
%                         value within some tolerance.
%   'DaviesBouldin'       The number of clusters giving the minimum
%                         DaviesBouldin value.
    OptimalK = nan;
       
%ClusteringFunction Clustering algorithm used.
%   The ClusteringFunction property is empty if the clustering solutions
%   are provided in the input. Otherwise it's a string or a function handle
%   representing the clustering algorithm requested.
    ClusteringFunction;
    
%NumObservations Number of observations.
%   The NumObservations property gives the number of observations used
%   when creating the object. This is the number of observations
%   supplied in the data, minus the number of rows with missing (NaN)
%   values.
    NumObservations;
    
     
%Missing A vector indicating which observations are missing.
%   The Missing property is a column vector of logical values, with one
%   value for each row in the input data. The value is TRUE if there are
%   missing (NaN) values in the corresponding row of the data matrix or
%   clustering solutions.
    Missing;
   
end


properties(GetAccess=protected, SetAccess=protected) 
   FunLoc=[];
    % An integer saving the location of builtin clustering algorithm 
    % in  {'kmeans', 'linkage', 'gmdistribution'} if Y is a string;
    % Otherwise it's empty.
    % 1 - represents "kmeans". 
    % 2 - represents "linkage"
    % 3 - represents "gmdistribution"
           
   ClustFuncHandle = [];
   % Either a function handle passed from the input, or a function handle 
   % based on the string representing the clustering algorithm
   % ClusteFuncHandle will be empty if "linkage" clustering option is
   % requested or a pre-computed clustering solution is provided.
   
   LinkageFuncHandle = []; 
   % If "linkage" clustering option is requested. LinkageFuncHandle will be
   % "linkage"; otherwise it will be empty
   
   ZToClustFuncHandle = []
   % If "linkage" clustering option is requested. ZToClustFuncHandle will be
   % "cluster"; otherwise it will be empty
   
   LinkageZ = [];
   % LinkageZ will save the Z output from LinkageFuncHandle if "linkage"
   % option is requested. Otherwise, LinkageZ will be empty.
   
   PrivX; 
    % This property contains the input X after removing the rows contains
    % NaN values.
    
   NewKList = [];
    %This property contains the new klist after removing the repeated
    %values appeared in both InspectedK property and new klist 
  
end

methods (Access = protected,Static) %(Hidden,Static)
    function  [args, dflts] =  getAllArgs()
        args  ={'klist'};
        dflts = {[]};
    end   
end

 methods (Hidden, Static)
       function a = empty(varargin)
            error(message('stats:clustering:evaluation:ClusterCriterion:NoEmptyAllowed'));
      end
 end
  
  
methods (Access = public)
   
    function this = addK(this,klist)
% ADDK Add additional number of clusters to be evaluated.
%    E = ADDK(E, KLIST) evaluates the number of clusters provided in KLIST
%    and then updates the related properties in clustering evaluation
%    object E. It can only be called when E.ClusteringFunction is not
%    empty. KLIST is a vector of positive integers containing the number of
%    clusters to be evaluated.
    
        if isempty(this.ClusteringFunction)
            error(message('stats:clustering:evaluation:ClusterCriterion:AddkforEmptyClusteringFunction'));
        end
        
        if isempty(this.X)
            error(message('stats:clustering:evaluation:ClusterCriterion:EmptyX'));
        end
        
        klist = checkKList(this,klist);
        
        % new check whether there is any overlapped between the old klist
        % and the new list. Remove the repeated K values if any
         unique_klist = unique([this.InspectedK klist],'stable');
         
         if numel(unique_klist) ~= numel([this.InspectedK klist])
               unique_klist(1:numel(this.InspectedK))=[];
               if isempty(unique_klist)
                    warning(message('stats:clustering:evaluation:ClusterCriterion:AllKHaveInspected'));
                    return;
               end
               warning(message('stats:clustering:evaluation:ClusterCriterion:RepeatedK'));
               this.NewKList = unique_klist;
         else
               this.NewKList = klist;
         end
   
      this.InspectedK =[this.InspectedK this.NewKList];
      
      if any(this.Missing)
          this.PrivX = this.X(~this.Missing,:);
      else
          this.PrivX = this.X;
      end
      
    end %addK
    
    
    function hout = plot(varargin)
%PLOT  Plot the requested criterion values 
%  PLOT(E) plots the requested criterion values versus the number of
%  clusters.
%
%  H = PLOT(E) returns the handle to the line in the plot.
%
%  H = PLOT(AX,E) plots into the axes AX.

        % Allow for PLOT(E), PLOT(AX,E), E.PLOT(AX)
        otype = 'clustering.evaluation.ClusterCriterion';
        [parent, this] = stats.internal.separateAxesAndObjectArguments(otype, varargin);
        hasParent = ~isempty(parent);
        narginchk(1,1+hasParent);

        [sortedKList, IDX] = sort(this.InspectedK);
        CritVals = this.CriterionValues;
        critVals = CritVals(IDX);

        [parent, ~, nextplot] = matlab.graphics.chart.internal.inputparsingutils.prepareAxes(parent, hasParent);
        h = plot(parent,sortedKList, critVals ,'-ob');
        ax  = ancestor(h,'axes');
        if ~isnan(this.OptimalK)
            optimalVal = this.CriterionValues(this.InspectedK == this.OptimalK);
            line(ax,this.OptimalK,optimalVal,'Marker','.','MarkerSize',10);
        end
        xlabel(ax,getString(message('stats:clustering:evaluation:ClusterCriterion:PlotXlabel')));
        YString = getString(message('stats:clustering:evaluation:ClusterCriterion:PlotYlabel',this.CriterionName));
        ylabel(ax,YString);
        parent.NextPlot = nextplot;
        if nargout>=1
            hout = h;
        end
    end
   
  function this = compact(this)
%COMPACT Compact this object.
%   C = COMPACT(E) sets E.X, E.OptimalY and E.Missing in clustering
%   evaluation object E to be empty.
  
      this.X = [];
      this.OptimalY = [];
      this.Missing = [];
  end
end

methods (Access = protected)
    function klist = checkKList(this,klist)
        %Sanity check for  klist
        if ~(isnumeric(klist) && isvector(klist) && all(klist > 0) ...
                && all(round(klist)==klist)) || any(isinf(klist)) || any(isnan(klist))
            error(message('stats:clustering:evaluation:ClusterCriterion:BadKList'));
        end
        
        badKs = find(klist>sum(~this.Missing,1));
        if ~isempty(badKs)
            klist(badKs) = [];
            if isempty(klist) 
                 error(message('stats:clustering:evaluation:ClusterCriterion:EmptyKLISTAfterRemoveBigKs'));
            end
            warning(message('stats:clustering:evaluation:ClusterCriterion:RemoveTooBigK'));
        end
        klist =klist(:)';
    end
function group = getPropertyGroups(obj)
 
  plist1 = struct('NumObservations',   obj.NumObservations,...
                  'InspectedK', obj.InspectedK,...
                  'CriterionValues',obj.CriterionValues,...
                  'OptimalK', obj.OptimalK);
  
  group(1) = matlab.mixin.util.PropertyGroup(plist1);

end

function footerText = getFooter(obj)
    % override used by the CustomDisplay interface to display the footer
    % use the same footer defined in
    % internal.stats.displayMethodsProperties for classreg objects
    footerText = [];
    internal.stats.displayMethodsProperties(obj);
end

function this = postProcess(this)
    %add NaNs back to the rows in OptimalY if there is missing values in
    %X or original clustering solution
    this.PrivX = [];
    if any(this.Missing) && ~isempty(this.OptimalY) && ...
            size(this.OptimalY,1) == this.NumObservations
        this.OptimalY = internal.stats.insertnan(this.Missing,this.OptimalY);
    end
end

function this = getKfromY(this, Y)
    %Y must be a matrix containing hard solutions
    nYCols = size(Y,2);
    for j = 1:nYCols
        this.InspectedK(j) = numel(unique(Y(~this.Missing,j)));
        %Warn if the number of clusters changes due to
        %removing rows with missing values
        oldK = sum(~isnan(unique(Y(:,j))));
        if oldK ~= this.InspectedK(j)
            warning(message('stats:clustering:evaluation:ClusterCriterion:KChangeInY',j,this.InspectedK(j)));
        end
    end
end

function IDX = soft2hard(this, softY)
%Convert cell array of matrix (soft solution) to matrix (hard solution)
    p = numel(this.InspectedK);
    IDX = zeros(size(this.X,1),p);
   
    for j = 1:p
        [~, IDX(:,j)] = max(softY{j},[],2);
        tempK = numel(unique(IDX(~this.Missing,j)));
        if tempK ~= this.InspectedK(j)
             warning(message('stats:clustering:evaluation:ClusterCriterion:Soft2hardBadNC',j,tempK));
             this.InspectedK(j) = tempK;
        end
        
    end
end

function IDX = evalFun(this,NC)
% Clustering PrivX into NC clustering using the requested clustering
% (either built-in clustering algorithm or a function handle)
       IDX = [];
       funcHandle = this.ClustFuncHandle;
       try
           if  ischar(this.ClusteringFunction) && this.FunLoc == 2% linkage clustering
               %this.LinkageZ must have already been computed before call this
               %method
               funcHandle = this.ZToClustFuncHandle;
               IDX = this.ZToClustFuncHandle(this.LinkageZ,NC);
            else
               IDX = this.ClustFuncHandle(this.PrivX,NC);
            end
       catch ME
           if strcmp('MATLAB:UndefinedFunction', ME.identifier) ...
                   && ~isempty(strfind(ME.message, func2str(funcHandle)))
               error(message('stats:clustering:evaluation:ClusterCriterion:FunNotFound', func2str(funcHandle)));
           else
               warning(message('stats:clustering:evaluation:ClusterCriterion:FunError', NC, func2str(funcHandle), ME.message));
               return;
           end
       end
       
       if isa(this.ClusteringFunction,'function_handle')
           %Sanity check on results from a function handle
           IDX = checkFuncHandle(this,IDX,NC,true);
       end
end

%Check whether the size if valid, convert a soft clustering solution
%into a hard clustering solution.
function funcOut = checkFuncHandle(this, funcOut,nc,isOriginal)
    if ~isnumeric(funcOut)
        error(message('stats:clustering:evaluation:ClusterCriterion:FuncOutNonNumeric'));
    end
    
    [nRows, nCols]=size(funcOut);
    if size(this.PrivX,1) ~= nRows
        error(message('stats:clustering:evaluation:ClusterCriterion:FuncOutSizeMismatch'));
    end
    if  any(any(isnan(funcOut)))
        error(message('stats:clustering:evaluation:ClusterCriterion:FuncOutNaN'));
    end
    if nc == 1
        if ~isvector(funcOut)
            error(message('stats:clustering:evaluation:ClusterCriterion:FuncOutBad'));
        end
        if numel(unique(funcOut)) == 1 && all(round(funcOut) == funcOut) && all(funcOut > 0)
            isHard = true;   %view funcOut as a hard clustering solution
        else
            isHard = false;  %view funcOut as a hard clustering solution even it may come from a soft clustering solution.
        end
        
    else % nc ~= 1
        if  nCols == 1  %vector
            if any(round(funcOut)~= funcOut) || any(funcOut <= 0)
                error(message('stats:clustering:evaluation:ClusterCriterion:FuncOutBad'));
            end
            actualNC =  numel(unique(funcOut));
            if actualNC ~= nc
                if isOriginal
                    warning(message('stats:clustering:evaluation:ClusterCriterion:FuncOutNCBadOnX',nc, actualNC));
                else
                    warning(message('stats:clustering:evaluation:ClusterCriterion:FuncOutNCBadOnRef',nc, actualNC));
                end
                funcOut = [];
                return;
            end
            isHard = true; %hard clustering solution
        else
            if nCols ~= nc
                error(message('stats:clustering:evaluation:ClusterCriterion:FuncOutBad'));
            end
            isHard = false; %soft score
        end
    end
    
    if ~isHard %convert soft solution to hard solution
        [~, funcOut] = max(funcOut,[],2);
        actualNC =  numel(unique(funcOut));
        if actualNC ~= nc
            if isOriginal
                warning(message('stats:clustering:evaluation:ClusterCriterion:SoftFuncOutNCBadOnX',nc,actualNC));
            else
                warning(message('stats:clustering:evaluation:ClusterCriterion:SoftFuncOutNCBadOnRef',nc, actualNC));
            end
            
            funcOut = [];
        end
    end
    
end

    %Base Class Constructor
    function this = ClusterCriterion(X,Y,klist)
               
         %Performs some basic sanity checking on X, Y, KList, also checks
         %missing values, finds the number of clusters to be evaluated, and
         %set the function handle for built-in algorithms
         
         if ~isfloat(X) || ~ismatrix(X)
             error(message('stats:clustering:evaluation:ClusterCriterion:XBad'));
         end
         if isempty(X)
             error(message('stats:clustering:evaluation:ClusterCriterion:XEmpty'));
         end
         
         this.X = X;
         this.Missing = any(isnan(this.X),2); %check missing values in X
       
         if isempty(klist)
             %Y must be precomputed clustering solutions if klist is empty
             if isnumeric(Y) && ismatrix(Y) %If Y contains hard clustering solution
                 [nYRows, nYCols] = size(Y);
                 if size(X,1) ~= nYRows;
                     error(message('stats:clustering:evaluation:ClusterCriterion:XHardYSizeMismatch'));
                 end
            
                 this.Missing = this.Missing | any(isnan(Y),2);
                 this.InspectedK = zeros(1,nYCols);
            
                 if ~(all(all(Y(~this.Missing,:) > 0)) && all(all(round(Y(~this.Missing,:))==Y(~this.Missing,:))))
                       error(message('stats:clustering:evaluation:ClusterCriterion:BadY'));
                 end
                 this = getKfromY(this, Y);

                                 
             elseif iscell(Y) %Soft clustering solution
                 p = numel(Y);
                 this.InspectedK=zeros(1,p);
                 for j=1:p
                     if ~ismatrix(Y{j}) || ~isnumeric(Y{j})
                          error(message('stats:clustering:evaluation:ClusterCriterion:BadY'));
                     end
                     [ny,K] = size(Y{j});
                     if ny ~= size(this.X,1)
                           error(message('stats:clustering:evaluation:ClusterCriterion:XSoftYSizeMismatch'));
                     end
                     this.Missing = this.Missing| any(isnan(Y{j}),2);
                     this.InspectedK(j) = K;      
                 end
                
             else 
                 error(message('stats:clustering:evaluation:ClusterCriterion:BadY'));
             end
             
             %We don't allow repeated K values if Y contains clustering
             %solutions.
             if numel(unique(this.InspectedK)) ~= numel(this.InspectedK)
                 error(message('stats:clustering:evaluation:ClusterCriterion:SameK'));
             end

         else  % If KLIST is provided, then Y must be a string or a function handle
             klist = checkKList(this,klist);
   
             unique_klist = unique(klist,'stable');
             if numel(unique_klist) ~= numel(klist) %Gives warning for repeated K values
                 warning(message('stats:clustering:evaluation:ClusterCriterion:RepeatedKinKList'));  
             end
             this.InspectedK = unique_klist;     
             
             if ischar(Y) %Y represents a built-in clustering algorithm
                 ClustFuncList = {'kmeans', 'linkage', 'gmdistribution'};
                 
                 [Y,this.FunLoc]= internal.stats.getParamVal(Y,ClustFuncList,'CLUST');
                 
                 if this.FunLoc == 1 %kmeans
                     this.ClustFuncHandle = ...
                         @(X,NC)(kmeans(X,NC, 'EmptyAction','singleton','rep',5));
                     
                 elseif this.FunLoc == 2
                     %Ward linkage with default Euclidean distance will be
                     %used. ClustFunctionHandle will not be used
                     this.LinkageFuncHandle =  @(X)(linkage(X,'ward'));
                     this.ZToClustFuncHandle = ...
                         @(Z,NC)(cluster(Z,'maxclust',NC));
                     
                 elseif this.FunLoc == 3
                     this.ClustFuncHandle = ...
                         @(X, NC)(cluster(fitgmdist(X,NC,'rep',5,'sharedcov',true),X));
                 end
                 
             elseif isa(Y, 'function_handle')
                 this.ClustFuncHandle = Y;
             else
                 error(message('stats:clustering:evaluation:ClusterCriterion:BadFun'));
             end
             
             this.ClusteringFunction = Y; %A string or a function handle
           
             
         end %if isempty(klist)
         
         this.NumObservations = sum(~this.Missing,1);
         
         if this.NumObservations == 0
             error(message('stats:clustering:evaluation:ClusterCriterion:NObservationsZero'));
         end
         
         if any(this.Missing)
             this.PrivX = this.X(~this.Missing,:);
         else
             this.PrivX = this.X;
         end
         
         this.CriterionValues = nan(1,numel(this.InspectedK));
    end
end % protected method


end% classdef


