 classdef (InferiorClasses = {?matlab.graphics.axis.Axes, ?matlab.ui.control.UIAxes}) GapEvaluation < clustering.evaluation.ClusterCriterion
%GapEvaluation Evaluation object based on the gap values
%   A GapEvaluation object computes the gap criterion values which can be
%   used for estimating the number of clusters. It's created using
%   the following syntax:
%
%   EVA = EVALCLUSTERS(X,CLUST,'GAP')
%
%   GapEvaluation properties:
%      X                      - Data used to create this object.
%      Distance               - The distance metric.
%      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 gap values for each number of clusters.
%      ClusteringFunction     - Clustering algorithm used.
%      Missing                - A vector indicating which observations are missing.
%      B                      - Number of reference data sets.
%      ReferenceDistribution  - Method used to generate reference data.
%      SearchMethod           - Method used to find the optimal clustering solution.
%      LogW                   - Natural logarithm of W on the input data.
%      ExpectedLogW           - Expectation of the natural logarithm of W on the reference data.
%      SE                     - Standard error of the expected LogW for each number of clusters on the reference data.
%      StdLogW                - Standard deviation of the natural logarithm of W on the reference data.
%
%   GapEvaluation methods:
%      addK             - Add a new list of number of clusters for evaluation.
%      compact          - Compact this object.
%      increaseB        - Add more simulation data sets.
%      plot             - Plot the gap values.
%
%   See also EVALCLUSTERS,clustering.evaluation.SilhouetteEvaluation
%            clustering.evaluation.CalinskiHarabaszEvaluation,
%            clustering.evaluation.DaviesBouldinEvaluation.

% Reference:
%    Tibshirani, R., Walther, G. and Hastie, T. (2001) Estimating the
%    number of data clusters via the Gap statistic. Journal of the Royal
%    Statistical Society B, 63, p411-423


%   Copyright 2013-2023 The MathWorks, Inc.
    
properties(GetAccess=public,SetAccess=protected)
%Distance The distance metric.
%   The Distance property specifies the distance metric that used. It's a
%   string representing a built-in distance measure, or a function handle.
    Distance = 'sqEuclidean';
    
%B The number of reference data sets.
%   The B property is an integer representing the number of reference data
%   sets used for computing gap criterion values.
    B = 100;
    
%ReferenceDistribution The method used to generate the reference data
%   The ReferenceDistribution property is a string indicating how reference
%   data are generated. It is 'uniform' or 'PCA'. The 'uniform' value means
%   generating each reference feature uniformly over the range of each
%   feature of the X property. 'PCA' means  generating the reference
%   feature from a uniform distribution over a box aligned with the
%   principal components of the X property.
    ReferenceDistribution = 'PCA'
       
%LogW Natural logarithm of W on the input data.
%   The LogW property is a numeric vector containing the natural logarithm
%   of W on the input data for each number of clusters in the InspectedK
%   property. W is the sum of the within-cluster average pairwise distance
%   values, summed over all the clusters.  If the distance is the squared
%   Euclidean distance, then W is also equal to the sum of within-cluster
%   sums of point-to-centroid squared Euclidean distances.
    LogW;
        
%ExpectedLogW Expectation of the natural logarithm of W on the reference data.
%   The ExpectedLogW property is a numeric vector containing the estimate
%   of expectation of the natural logarithm of W on the reference data for
%   each number of clusters in the InspectedK property. W is the sum of the
%   within-cluster average pairwise distances, summed over all the clusters.
%   If the distance is the squared Euclidean distance, then W is also the
%   sum of within-cluster sum of point-to-centroid squared Euclidean
%   distances.
    ExpectedLogW;
        
%StdLogW Standard deviation of the natural logarithm of W on the reference data.
%   The StdLogW property is a numeric vector containing the  standard
%   deviation of the natural logarithm of W on reference data for each
%   number of clusters. W is the sum of the within-cluster average pairwise
%   distance values, summed over all the clusters. If the distance is the
%   squared Euclidean distance, then W is also the sum of within-cluster
%   sums of point-to-centroid squared Euclidean distances.
    StdLogW;
        
%SE Standard error of the natural logarithm of W on the reference data.
%   The SE property is a numeric vector containing the standard error of
%   the natural logarithm of W on reference data for each number of
%   clusters in the InspectedK property. W is the sum of the within-cluster
%   average pairwise distance values, summed over all the clusters. If the
%   distance is the squared Euclidean distance, then W is also the sum of
%   within-cluster sum of point-to-centroid squared Euclidean distances.
    SE; 
        
%SearchMethod Method used to find the optimal clustering solution.
%   The SearchMethod property is a string indicating how the optimal number
%   of clusters is estimated. It's either 'FirstMaxSe' or 'GlobalMaxSe'.
    SearchMethod; 
 end
    
    properties(GetAccess=protected,SetAccess=protected)
        VTransform % the transfer of V, [U,S,V]=svd(X0)
        minX2;  % The minimal value of X2, X2= X0*V; it will be empty if PCA is used
        maxX2;  % The maximal value of X2, X2= X0*V; it will be empty if PCA is used
        refDistLoc; %1:PCA reference distribution, 2: uniform reference distribution
        SearchMethodLoc; %1-FirstMaxSe 2 -GlobalMaxSE
        NXCols; % The number of columns in X
        minX;
        maxX;
        needSort = false; %True means the input klist are not sorted
        distanceLoc;
        distLoc;
        ValidResult;
    end
    
    methods
        function this = addK(this,klist)
            
            this = addK@clustering.evaluation.ClusterCriterion(this,klist);
            klist = this.NewKList;
            NumNewK = numel(klist);
            logW = nan(1,NumNewK);
            criterionValues = nan(1,NumNewK);
            goodResult = true(1,NumNewK);
            [refLogW,goodResult] = getRefLogW(this,this.B,klist,goodResult);
            [expLogW,stdLogWRef,SERef] = getRefStats(this,refLogW);
                   
            NumOldK = numel(this.CriterionValues);
            isNewK = [false(1,NumOldK) true(1,NumNewK)];
            
            this.CriterionValues = [this.CriterionValues, criterionValues];
            this.StdLogW =[this.StdLogW, stdLogWRef];
            this.ExpectedLogW = [this.ExpectedLogW, expLogW];
            this.SE = [this.SE, SERef];
            this.LogW = [this.LogW logW];
            this.ValidResult = [this.ValidResult goodResult];
            %Sort the K values if needed
            [sortedKlist,sortIdx] = sort(this.InspectedK);
            
            if this.needSort || any(sortedKlist ~= this.InspectedK)
                this.needSort = true;
                this = sortLists(this,sortedKlist,sortIdx);
                isNewK = isNewK(sortIdx); 
            end
            
            %Clustering input data using the K values in the new KList,
            p = numel(this.InspectedK);
            oldOptimalK = this.OptimalK;
            oldOptimalY = this.OptimalY;
            this.OptimalK = NaN;
            this.OptimalY = [];
            if this.SearchMethodLoc == 1 %FirstMaxSE
                for j = p:-1:1
                    NC = this.InspectedK(j);
                    IDX = [];
                    if isNewK(j) 
                        %If this K value is new, evaluate this value
                        [this,IDX] = getGapValue(this,j);
                    end
                    if ~this.ValidResult(j)
                           continue;
                    end
                    %find the next valid clustering result
                    nextValid = find(this.ValidResult(j+1:end),1);
                   
                    if  isempty(nextValid) ... %no next valid clustering result
                           || (this.CriterionValues(j) > 0 && this.CriterionValues(j) >=...
                                (this.CriterionValues(j+nextValid)-this.SE(j+nextValid)))       
                        this = updateOptimal(this,isNewK(j),NC,oldOptimalY,IDX,NC==oldOptimalK);
                    end
                end
                
            else%GlobalMaxSE
                maxGap= -inf;
                th = -inf;
                       
                for j=p:-1:1
                    NC = this.InspectedK(j);
                    IDX = [];
                    if isNewK(j) %this is a new K value
                        [this,IDX] = getGapValue(this,j);
                    end
                    if this.CriterionValues(j) >= maxGap
                        maxGap = this.CriterionValues(j);
                        th = maxGap - this.SE(j);
                    end
                    if  this.CriterionValues(j) >= th && this.CriterionValues(j) > 0
                        this = updateOptimal(this,isNewK(j),NC,oldOptimalY,IDX,NC==oldOptimalK);
                    end
                end
            end
            
            %Fill the OptimalY property if it's empty
            if ~isnan(this.OptimalK) && isempty(this.OptimalY)
                if ~(ischar(this.ClusteringFunction) && this.FunLoc == 2)
                    warning(message('stats:clustering:evaluation:GapEvaluation:NonOriginalOptimalY'));
                end
                this.OptimalY = this.evalFun(this.OptimalK);
           end
            
            %Re-sort some properties based on the original order of klist
            this = backToOriginalList(this,sortIdx);
            this = postProcess(this);
            this.NewKList = [];
           
      end % addK
      
      function this = increaseB(this, B2)
%INCREASEB increase the number of reference data sets
%   E = INCREASEB(E, B) generates B more reference data sets and updates
%   the gap criterion values and other related properties in the
%   GapEvaluation object E.
          
          checkB(B2); %Sanity check on B2
          if isempty(this.X)
            error(message('stats:clustering:evaluation:ClusterCriterion:EmptyX'));
          end
          if any(this.Missing)
              this.PrivX = this.X(~this.Missing,:);
          else
              this.PrivX = this.X;
          end
          
          [sortedKlist,sortIdx]= sort(this.InspectedK);
          if this.needSort
              this = sortLists(this,sortedKlist,sortIdx);
          end
          this.B = this.B+B2;
          %Clustering B simulation data sets and compute all the statistics
                 
          [refLogW, this.ValidResult] = getRefLogW(this,B2,this.InspectedK, this.ValidResult);
          this = IncreUpdateRefStats(this,refLogW);
          %Clustering input data set and find the optimal solution.
          this = findOptimal(this);
          %Resort klist and CriterionValues based on the original order
          this = backToOriginalList(this,sortIdx);
          
          this = postProcess(this);
      end % increaseB
      
      function hout = plot(varargin)
%PLOT  Plot the gap values.
%  PLOT(E) plots the gap values with one standard error bars versus the
%  number of clusters.
%
%  H = PLOT(E) returns the handle to the line in the plot.

            otype = 'clustering.evaluation.ClusterCriterion';
            [parent, this] = stats.internal.separateAxesAndObjectArguments(otype, varargin);
            hasParent = ~isempty(parent);
            narginchk(1,1+hasParent);
            [parent, ~, nextplot] = matlab.graphics.chart.internal.inputparsingutils.prepareAxes(parent, hasParent);
            
            if this.needSort
                [sortedKList, IDX] = sort(this.InspectedK);
                h = errorbar(parent,sortedKList,this.CriterionValues(IDX),this.SE(IDX),'-ob');
            else
                h = errorbar(parent,this.InspectedK,this.CriterionValues,this.SE,'-ob');
            end
            
            if ~isnan(this.OptimalK)
                optimalVal = this.CriterionValues(this.InspectedK == this.OptimalK);
                line(parent,this.OptimalK,optimalVal,'Marker','.','MarkerSize',10);
            end
            xlabel(parent,getString(message('stats:clustering:evaluation:ClusterCriterion:PlotXlabel')));
            YString = getString(message('stats:clustering:evaluation:ClusterCriterion:PlotYlabel',this.CriterionName));
            ylabel(parent,YString);
            parent.NextPlot = nextplot;
            if nargout > 0
                hout = h;
            end
        end
    end
    
   methods (Hidden)
        function   this = GapEvaluation(X,Y,varargin)
            narginchk(2,inf);
            [basePnames, baseDflts ]= clustering.evaluation.ClusterCriterion.getAllArgs();
            
            [klist,~,args] = ...
                internal.stats.parseArgs(basePnames, baseDflts, varargin{:});
                           
            if isempty(klist) %klist must be provided for the gap criterion
                error(message('stats:clustering:evaluation:GapEvaluation:EmptyKList'));
            end
            
            this = this@clustering.evaluation.ClusterCriterion(X,Y,klist);
            this.CriterionName = 'Gap';
               
            defaultSearch ='GlobalMaxSE';

            pnames = {'Distance'     'ReferenceDistribution'  'B' 'SearchMethod'};
            dflts =  {'sqeuclidean'  'PCA'                    100  defaultSearch };
            
            [distance, refDist,BNum, searchMethod]=...
                internal.stats.parseArgs(pnames, dflts, args{:});
            
            %Sanity check for distance
            %Distance must be a string or a function handle
            distList = {'sqEuclidean' 'Euclidean' 'cityblock'  'cosine' 'correlation'};
            if ischar(distance)
                [this.Distance,this.distLoc] = internal.stats.getParamVal(...
                    distance,distList,'distance');
            elseif isa(distance,'function_handle')
                this.distLoc = 0;
                this.Distance = distance;
            else
                error(message('stats:clustering:evaluation:GapEvaluation:BadDist'));
            end
           
            if this.distLoc ~= 1 && ischar(Y)
                %if Y is a built-in function and the distance is not the sqEuclidean
                %distance
                if this.FunLoc == 1
                    this.ClustFuncHandle = getKmeansFunc(this.Distance);
                    %The above line is the same as the following code: 
                    %this.ClustFuncHandle = ...
                    %   @(X,NC)(kmeans(X,NC,'rep',5,'EmptyAction','singleton','distance',this.Distance));
                    %We assign the function handle inside the sub-function
                    %to make sure the onCleanup object destructor will be
                    %called when needed.
                    
                elseif this.FunLoc == 2 %linkage
                    if this.distLoc ~= 2 %the distance is neither 'euclidean' nor 'sqeuclidean'
                        this.LinkageFuncHandle = getLinkageFuncHandle(this.Distance);
                        %The above line is same as the following code.
                        %We assign the function handle inside the sub-function
                        %to make sure the onCleanup object destructor will be
                        %called when needed.
                        % this.LinkageFuncHandle = ...
                        %   @(X)(linkage(X,'linkage','ave',this.Distance);
                     end
                   
                elseif this.FunLoc == 3 %gmdistribution
                    error(message('stats:clustering:evaluation:GapEvaluation:InvalidDistForGM'));
                end
            end
            
            %Sanity check for ReferenceDistribution
            [this.ReferenceDistribution,this.refDistLoc]=internal.stats.getParamVal(...
                refDist,{'PCA','Uniform'},'ReferenceDistribution');
            
            if this.refDistLoc == 1; %PCA reference distribution
                this = getSVD(this);
            else
                this.maxX = max(this.PrivX,[],1);
                this.minX = min(this.PrivX,[],1);
            end
            
            checkB(BNum);
            this.B = BNum;
            
            [this.SearchMethod,this.SearchMethodLoc] = ...
                internal.stats.getParamVal(...
                searchMethod,{'FirstMaxSE','GlobalMaxSE'},'SearchMethod'); 
            
            this.NXCols = size(X,2);
            
           [sortedKlist,sortIdx] = sort(this.InspectedK);
            
            if any(sortedKlist ~= this.InspectedK)
                this.needSort = true;
                this.InspectedK = sortedKlist;
            end
            
            %Initialize the properties:
            p = numel(this.InspectedK);
            this.LogW = nan(1,p);
            this.CriterionValues = nan(1,p);
            this.ExpectedLogW = nan(1,p);
            this.StdLogW = nan(1,p);
            this.SE = nan(1,p);
            this.ValidResult = true(1,p);
            
            %Clustering B simulation data sets and compute all the
            %statistics on the reference data sets.
            [refLogW,this.ValidResult] = getRefLogW(this,this.B,this.InspectedK, this.ValidResult ); 
            [this.ExpectedLogW,this.StdLogW,this.SE] = getRefStats(this,refLogW);
            %Pre-compute the Z matrix for input data if the linkage option is requested
            if ~isempty(this.FunLoc) && this.FunLoc == 2
                    this.LinkageZ = this.LinkageFuncHandle(this.PrivX);   
            end
               
            %Clustering input data after the clustering on reference data
            %is done. Find the optimal number of clusters.
            this = findOptimal(this);
                        
            %Re-sort some properties based on the original order of input klist
            this = backToOriginalList(this,sortIdx);
            
            this = postProcess(this);
        end %Constructor
        
   end 
    
    
   methods (Access = protected)
       function IDX = evalFunXref(this,NC,Xref)
           %clustering the reference data using the original reference data
           try
               IDX = this.ClustFuncHandle(Xref,NC);
           catch ME
               if strcmp('MATLAB:UndefinedFunction', ME.identifier) ...
                       && ~isempty(strfind(ME.message, func2str(this.ClustFuncHandle)))
                   error(message('stats:clustering:evaluation:ClusterCriterion:FunNotFound', func2str(this.ClustFuncHandle )));
               else
                   warning(message('stats:clustering:evaluation:GapEvaluation:FunErrorOnRef',NC, func2str(this.ClustFuncHandle), ME.message));
                   IDX = [];
                   return;
               end
           end
           if isa(this.ClusteringFunction,'function_handle')
               IDX = checkFuncHandle(this,IDX,NC,false);
           end
           
       end
        
        function IDX = evalFunZref(this,NC,Zref)
            %clustering the reference data using the Z matrix generated from the original reference data
            try
                IDX = this.ZToClustFuncHandle(Zref,NC);
            catch ME
                if strcmp('MATLAB:UndefinedFunction', ME.identifier) ...
                        && ~isempty(strfind(ME.message, func2str(this.ZToClustFuncHandle)))
                    error(message('stats:clustering:evaluation:ClusterCriterion:FunNotFound', func2str(this.ZToClustFuncHandle)));
                else
                    warning(message('stats:clustering:evaluation:GapEvaluation:FunErrorOnRef',NC, func2str(this.ZToClustFuncHandle), ME.message));
                    IDX = [];
                end
            end
        end
        
        function this = findOptimal(this)
            %Run clustering algorithm with the number specified in
            %this.InspectedK. Update the corresponding properties and
            %update the OptimalY and the OptimalK properties.
            p = numel(this.InspectedK);
            this.OptimalK = NaN;
            this.OptimalY = [];
            if this.SearchMethodLoc == 1 %FirstMaxSE
                   for j = p:-1:1
                       NC = this.InspectedK(j);
                       [this,IDX] = getGapValue(this,j);
                      
                       if ~this.ValidResult(j)
                           continue;
                       end
                       nextValid = find(this.ValidResult(j+1:end),1);
                       if  isempty(nextValid) ...
                           || (this.CriterionValues(j) > 0 && this.CriterionValues(j) >=...
                                (this.CriterionValues(j+nextValid)-this.SE(j+nextValid)))
                            this.OptimalK = NC;
                            this.OptimalY = IDX;
                       end
                        
                    end
                   
                else%GlobalMaxSE
                    maxGap= -inf;
                    th = -inf;
            
                    for j=p:-1:1
                        NC = this.InspectedK(j);
                        [this,IDX] = getGapValue(this,j);
                        if this.CriterionValues(j) >= maxGap
                            maxGap = this.CriterionValues(j);
                            th = maxGap - this.SE(j);
                        end
                        if  this.CriterionValues(j) >= th && this.CriterionValues(j) > 0
                            this.OptimalK = NC;
                            this.OptimalY = IDX;
                        end        
                    end
            end
        end
        
        function [this,IDX] = getGapValue(this,j)
            %Update the gap value and LogW value for the jth value in the
            %sorted this.InspectedK 
              if ~this.ValidResult(j)
                %no need to clustering X with this particular number of
                %clusters because clustering on the reference data has been
                %failed at least once.
                 IDX = [];
                 this.CriterionValues(j) = nan;
                 this.LogW(j) = nan;
                 return;
              end
            
            NC = this.InspectedK(j);
            if ischar(this.ClusteringFunction) && this.FunLoc==1 && this.distLoc==1
                [IDX, ~, SUMD] = kmeans(this.PrivX,NC,'rep',5,'empty','singleton');
                this.LogW(j) = log(sum(SUMD,1));
            else
                IDX = this.evalFun(NC);
                if ~isempty(IDX)
                    this.LogW(j) = getLogW(this.PrivX,IDX,this.Distance,this.distLoc==1);
                else
                    this.LogW(j)= nan;
                    this.ValidResult(j) = false;
                end
            end
            
            %Compute the gap value
            this.CriterionValues(j) = this.ExpectedLogW(j)-this.LogW(j);
        end
        
        function [ExpLogW,StdLogW,SE] = getRefStats(this,refLogW)
            %Get the Expected value,STD, SE and the variance of logW after
            %refLogW are computed
            %Any columns that has NaN values will make the corresponding
            %entries in the following row vector as NaN.
            ExpLogW = mean(refLogW,1);
            StdLogW = sqrt(var(refLogW,1,1)); 
            SE = StdLogW*sqrt(1/this.B+1); %The standard error
        end
          
        function [refLogW, validResult] = getRefLogW(this,B,klist, validResult)
            %validResult is a vector recording which entries are valid
            if ischar(this.ClusteringFunction) && this.FunLoc==1 && this.distLoc == 1
                 warnCell={ 'stats:kmeans:EmptyClusterInBatchUpdate',...
                'stats:kmeans:FailedToConvergeRep',...
                'stats:kmeans:EmptyClusterRep'};
            else
                 warnCell={'stats:gmdistribution:FailedToConverge',...
                'stats:gmdistribution:FailedToConvergeReps',...
                'stats:gmdistribution:IllCondCov',...
                'stats:kmeans:EmptyClusterInBatchUpdate',...
                'stats:kmeans:EmptyCluster',...
                'stats:kmeans:EmptyClusterRep',...
                'stats:kmeans:FailedToConverge',...
                'stats:kmeans:FailedToConvergeRep'};
            end
            
            p = numel(klist);
            refLogW = nan(B,p);
            
            for whichID = 1:numel(warnCell)
                status(whichID) = warning('off', warnCell{whichID});
            end
            restoreWS = onCleanup(@() warning(status));
            
            for b = 1:B
                %Generate reference data B times
                XRef = getReferenceData(this);
                ZRef = [];
                if ~isempty(this.FunLoc) && this.FunLoc == 2
                    ZRef = this.LinkageFuncHandle(XRef);
                end
                for j = 1:p
                    if ~validResult(j)
                         %If clustering reference data into NC clusters has
                         %failed before, there is no need to clustering on other
                         %reference data sets with the same NC
                         continue;      
                    end
                    NC = klist(j);
                    if ischar(this.ClusteringFunction) && this.FunLoc==1 && this.distLoc == 1
                        [~, ~,SUMD_ref] = kmeans(XRef,NC,'rep',5,'empty','singleton');
                        refLogW(b,j) = log(sum(SUMD_ref,1));
                    else% other clustering except kmeans
                        if isempty(ZRef)
                            IDX = this.evalFunXref(NC,XRef);
                        else
                            IDX = this.evalFunZref(NC,ZRef);
                        end
                        if ~isempty(IDX)
                            refLogW(b,j) = getLogW(XRef,IDX,this.Distance,this.distLoc==1);
                        else
                            validResult(j) = false;  % refLogW(b,j) will remain as NaN values
                        end
                    end
                end
            end
            warning(status);
        end
        

        function  this = IncreUpdateRefStats(this,newRefLogW)
            %Incrementally update the ExpectedLogW the StdLogW
            %and the SE property after getting the logW matrix on new reference data 
             B2 = size(newRefLogW,1);
             newSumLogW = sum(newRefLogW,1);
             newSSELogW = var(newRefLogW,1,1)*B2; 
             B1 = this.B-B2;
             oldSumLogW = this.ExpectedLogW*B1;
             oldSSELogW = (this.StdLogW.^2)*B1;
             SSELogW = oldSSELogW+newSSELogW+(B1/((this.B)*B2))*(((B2/B1)*oldSumLogW-newSumLogW).^2);
             this.ExpectedLogW = (newSumLogW+ oldSumLogW)/(this.B);
             this.StdLogW = sqrt(SSELogW/this.B);
             this.SE = this.StdLogW*sqrt(1/this.B+1);
       
        end
        
                
        %Generate reference data
        function Xref = getReferenceData(this)
            if this.refDistLoc == 1%PCA reference
                Xref = getUniformRef(this.minX2,this.maxX2,...
                    this.NumObservations,this.NXCols);
                Xref = Xref*this.VTransform;
                %The mean of Xref has been shifted, but it's OK because
                %shifting the mean doesn't affect the separation of the
                %data.  
            else % uniform reference
                Xref = getUniformRef(this.minX,this.maxX, this.NumObservations,...
                    this.NXCols);
            end
        end
        
        function this = getSVD(this)
            meanX=mean(this.PrivX);
            X0= bsxfun(@minus,this.PrivX,meanX);
            [~,~,V] = svd(X0);
            X2=X0*V;
            this.minX2 = min(X2,[],1);
            this.maxX2 = max(X2,[],1);
            this.VTransform = V';
        end
       
        function this = backToOriginalList(this,sortIdx)
            if this.needSort
                this.InspectedK(sortIdx) = this.InspectedK;
                this.CriterionValues(sortIdx) = this.CriterionValues;
                this.LogW(sortIdx) = this.LogW;
                this.StdLogW(sortIdx) = this.StdLogW;
                this.ExpectedLogW(sortIdx) = this.ExpectedLogW;
                this.SE(sortIdx) = this.SE;
                this.ValidResult(sortIdx) = this.ValidResult;
            end
            this.CriterionValues(~this.ValidResult) = nan;
            this.LogW(~this.ValidResult) = nan;
            this.StdLogW(~this.ValidResult) = nan;
            this.ExpectedLogW(~this.ValidResult) = nan;
            this.SE(~this.ValidResult) = nan;        
            
        end
        
        % Sort some properties based on sortIdx
        function this = sortLists(this,sortedKlist,sortIdx)
            this.InspectedK = sortedKlist;
            this.CriterionValues =this.CriterionValues(sortIdx);
            this.StdLogW = this.StdLogW(sortIdx);
            this.ExpectedLogW = this.ExpectedLogW(sortIdx);
            this.SE = this.SE(sortIdx);
            this.LogW = this.LogW(sortIdx);
            this.ValidResult = this.ValidResult(sortIdx);
        end
        
        %Update OptimalK and OptimalY for addK method
        function this = updateOptimal(this,isNewK,NC,oldOptimalY,IDX,isOldOptimalK)
            this.OptimalK = NC;
            this.OptimalY = IDX;
            if ~isNewK && isOldOptimalK
                 this.OptimalY = oldOptimalY;
            end
        end
    end%protected methods
    
  end%class


function  logW = getLogW(X,index,dist,isSqEuclidean)
% Get log(W), where W is explained in formula 2 in Reference 
p = size(X,2);
clusts = unique(index); 
num = numel(clusts);
centroids = nan(num,p);
sumD = nan(num,1);

if isSqEuclidean
    %get sum of squared Euclidean distance of point-to-centroid distances
    %if the distance is the squared Euclidean distance
    for i = 1:num
        members = (index == clusts(i));
        if any(members)
            ni = sum(members);
            centroids(i,:) = sum(X(members,:),1) / ni;
            sumD(i)= sum(pdist2(X(members,:),centroids(i,:)).^2);
            
        end
    end
else %other distance rather than squared Euclidean distance
    %Get the sum of pairwise distances divided the number of points 
    for i = 1:num
        members = (index == clusts(i));
        if any(members)
            ni = sum(members);
            sumD(i) = sum(pdist(X(members,:),dist))/ni;
        end
    end
end

logW = log(sum(sumD));
end

%Generate uniform data between on the interval from minX and maxX
function Xref = getUniformRef(minX,maxX,nXRows,nXCols)
  Xref = bsxfun(@plus,minX,bsxfun(@times,(maxX-minX),rand(nXRows,nXCols)));
end

%Sanity check on the number of simulations 
function checkB(B)
     if ~isscalar(B)||~isnumeric(B) || B <= 0 || round(B) ~= B || isinf(B) ||  isnan(B)
          error(message('stats:clustering:evaluation:GapEvaluation:InvalidB'));    
     end
end

function fh = getKmeansFunc(distMetric)
%fh = @(X,NC)(kmeans(X,NC,'rep',5,'EmptyAction','singleton','distance',this.Distance));
     fh = @nested;
     function out = nested(X, NC)
        out=kmeans(X,NC, 'rep',5,'EmptyAction','singleton','distance',distMetric);
      end
end


function fh = getLinkageFuncHandle(distMetric)
% fh = @(X,NC)(clusterdata(X,'linkage','ave',this.Distance));
     fh = @nested;
     function out = nested(X)
        out= linkage(X,'ave',distMetric);
      end
end