# HG changeset patch
# User Stefano Guidoni <ilguido@users.sf.net>
# Date 1617465072 -7200
#      Sat Apr 03 17:51:12 2021 +0200
# Node ID 933bc06e11e465f6c29fcf40e9386dbb69d7211a
# Parent  d354b1b49cad43c73fdcfdf53034d137a91f3ffa
evalclusters: new function

diff -r d354b1b49cad -r 933bc06e11e4 INDEX.in
--- a/INDEX.in	Sun Mar 28 19:30:40 2021 +0200
+++ b/INDEX.in	Sat Apr 03 17:51:12 2021 +0200
@@ -104,6 +104,7 @@
 Clustering
  cmdscale
  cophenet
+ evalclusters
  kmeans
  linkage
  mahal
diff -r d354b1b49cad -r 933bc06e11e4 NEWS
--- a/NEWS	Sun Mar 28 19:30:40 2021 +0200
+++ b/NEWS	Sat Apr 03 17:51:12 2021 +0200
@@ -16,6 +16,8 @@
  
  ** datasample         (patch #10046)
  
+ ** evalclusters
+ 
  ** gscatter           (patch #10043)
 
  ** inconsistent       (patch #10008)
diff -r d354b1b49cad -r 933bc06e11e4 inst/evalclusters.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/evalclusters.m	Sat Apr 03 17:51:12 2021 +0200
@@ -0,0 +1,352 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{eva} =} evalclusters (@var{x}, @var{clust}, @var{criterion})
+## @deftypefnx {Function File} {@var{eva} =} evalclusters (@dots{}, @qcode{Name}, @qcode{Value})
+##
+## Create a clustering evaluation object to find the optimal number of clusters.
+##
+## @code{evalclusters} creates a clustering evaluation object to evaluate the
+## optimal number of clusters for data @var{x}, using criterion @var{criterion}.
+## The input data @var{x} is a matrix with @code{n} observations of @code{p}
+## variables.
+## The evaluation criterion @var{criterion} is one of the following:
+## @table @code
+## @item @qcode{CalinskiHarabasz} 
+## to create a @code{CalinskiHarabaszEvaluation} object.
+##
+## @item @qcode{DaviesBouldin} 
+## to create a @code{DaviesBouldinEvaluation} object.
+##
+## @item @qcode{gap} 
+## to create a @code{GapEvaluation} object.
+##
+## @item @qcode{silhouette} 
+## to create a @code{SilhouetteEvaluation} object.
+##
+## @end table
+## The clustering algorithm @var{clust} is one of the following:
+## @table @code
+## @item @qcode{kmeans} 
+## to cluster the data using @code{kmeans} with @code{EmptyAction} set to
+## @code{singleton} and @code{Replicates} set to 5.
+##
+## @item @qcode{linkage} 
+## to cluster the data using @code{clusterdata} with @code{linkage} set to
+## @code{Ward}.
+##
+## @item @qcode{gmdistribution} 
+## to cluster the data using @code{fitgmdist} with @code{SharedCov} set to
+## @code{true} and @code{Replicates} set to 5.
+##
+## @end table
+## If the @var{criterion} is @code{CalinskiHarabasz}, @code{DaviesBouldin}, or
+## @code{silhouette}, @var{clust} can also be a function handle to a function
+## of the form @code{c = clust(x, k)}, where @var{x} is the input data,
+## @var{k} the number of clusters to evaluate and @var{c} the clustering result.
+## The clustering result can be either an array of size @code{n} with @code{k}
+## different integer values, or a matrix of size @code{n} by @code{k} with a
+## likelihood value assigned to each one of the @code{n} observations for each
+## one of the @var{k} clusters. In the latter case, each observation is assigned
+## to the cluster with the higher value.
+## If the @var{criterion} is @code{CalinskiHarabasz}, @code{DaviesBouldin}, or
+## @code{silhouette}, @var{clust} can also be a matrix of size @code{n} by 
+## @code{k}, where @code{k} is the number of proposed clustering solutions, so
+## that each column of @var{clust} is a clustering solution.
+##
+## In addition to the obligatory @var{x}, @var{clust} and @var{criterion} inputs
+## there is a number of optional arguments, specified as pairs of @code{Name} 
+## and @code{Value} options. The known @code{Name} arguments are:
+## @table @code
+## @item @qcode{KList} 
+## a vector of positive integer numbers, that is the cluster sizes to evaluate.
+## This option is necessary, unless @var{clust} is a matrix of proposed 
+## clustering solutions.
+##
+## @item @qcode{Distance} 
+## a distance metric as accepted by the chosen @var{clust}. It can be the
+## name of the distance metric as a string or a function handle. When
+## @var{criterion} is @code{silhouette}, it can be a vector as created by 
+## function @code{pdist}. Valid distance metric strings are: @code{sqEuclidean}
+## (default), @code{Euclidean}, @code{cityblock}, @code{cosine}, 
+## @code{correlation}, @code{Hamming}, @code{Jaccard}.
+## Only used by @code{silhouette} and @code{gap} evaluation.
+##
+## @item @qcode{ClusterPriors} 
+## the prior probabilities of each cluster, which can be either @code{empirical}
+## (default), or @code{equal}. When @code{empirical} the silhouette value is
+## the average of the silhouette values of all points; when @code{equal} the
+## silhouette value is the average of the average silhouette value of each
+## cluster. Only used by @code{silhouette} evaluation.
+##
+## @item @qcode{B} 
+## the number of reference datasets generated from the reference distribution.
+## Only used by @code{gap} evaluation.
+##
+## @item @qcode{ReferenceDistribution} 
+## the reference distribution used to create the reference data. It can be
+## @code{PCA} (default) for a distribution based on the principal components of
+## @var{X}, or @code{uniform} for a uniform distribution based on the range of
+## the observed data. @code{PCA} is currently not implemented.
+## Only used by @code{gap} evaluation.
+##
+## @item @qcode{SearchMethod} 
+## the method for selecting the optimal value with a @code{gap} evaluation. It 
+## can be either @code{globalMaxSE} (default) for selecting the smallest number
+## of clusters which is inside the standard error of the maximum gap value, or
+## @code{firstMaxSE} for selecting the first number of clusters which is inside
+## the standard error of the following cluster number.
+## Only used by @code{gap} evaluation.
+##
+## @end table
+##
+## Output @var{eva} is a clustering evaluation object.
+##
+## @end deftypefn
+##
+## @seealso{CalinskiHarabaszEvaluation, DaviesBouldinEvaluation, GapEvaluation, 
+## SilhouetteEvaluation}
+
+function cc = evalclusters (x, clust, criterion, varargin)
+  ## input check
+  if (nargin < 3)
+    print_usage ();
+  endif
+  
+  ## parsing input data
+  if ((! ismatrix (x)) || (! isnumeric (x)))
+    error ("evalclusters: 'x' must be a numeric matrix");
+  endif
+  
+  ## useful values for input check
+  n = rows (x);
+  p = columns (x);
+  
+  ## parsing the clustering algorithm
+  if (ischar (clust))
+    clust = lower (clust);
+    if (! any (strcmp (clust, {"kmeans", "linkage", "gmdistribution"})))
+      error ("evalclusters: unknown clustering algorithm '%s'", clust);
+    endif
+  elseif (! isscalar (clust))
+    if ((! isnumeric (clust))  || (length (size (clust)) != 2) || ...
+        (rows (clust) != n))
+      error ("evalclusters: invalid matrix of clustering solutions");
+    endif
+  elseif (! isa (clust, "function_handle"))
+    error ("evalclusters: invalid argument for 'clust'");
+  endif
+  
+  ## parsing the criterion parameter
+  ## we check the rest later, as the check depends on the chosen criterion
+  if (! ischar (criterion))
+    error ("evalclusters: invalid criterion, it must be a string");
+  else
+    criterion = lower (criterion);
+    if (! any (strcmp (criterion, {"calinskiharabasz", "daviesbouldin", ...
+                                   "silhouette", "gap"})))
+      error ("evalclusters: unknown criterion '%s'", criterion);
+    endif
+  endif
+  
+  ## some default value
+  klist = [];
+  distance = "sqeuclidean";
+  clusterpriors = "empirical";
+  b = 100;
+  referencedistribution = "pca";
+  searchmethod = "globalmaxse";
+  
+  ## parse the name/value pairs
+  pair_index = 1;
+  while (pair_index < (nargin - 3))
+    ## type check
+    if (! ischar (varargin{pair_index}))
+      error ("evalclusters: invalid property, string expected");
+    endif
+    
+    ## now parse the parameter
+    switch (lower (varargin{pair_index}))
+      case "klist"
+        ## klist must be an array of positive interger numbers;
+        ## there is a special case when it can be empty, but that is not the
+        ## suggested way to use it (it is better to omit it instead)
+        if (isempty (varargin{pair_index + 1}))
+          if (ischar (clust) || isa (clust, "function_handle"))
+            error (["evalclusters: 'KList' can be empty only when 'clust' "...
+                    "is a matrix"]);
+          endif
+        elseif ((! isnumeric (varargin{pair_index + 1})) || ...
+            (! isvector (varargin{pair_index + 1})) || ...
+            any (find (varargin{pair_index + 1} < 1)) || ...
+            any (floor (varargin{pair_index + 1}) != varargin{pair_index + 1}))
+          error ("evalclusters: 'KList' must be an array of positive integers");
+        endif
+        klist = varargin{pair_index + 1};
+      case "distance"
+        ## used by silhouette and gap
+        if (! (strcmp (criterion, "silhouette") || strcmp (criterion, "gap")))
+          error (["evalclusters: distance metric cannot be used with '%s'"...
+                  " criterion"], criterion);
+        endif
+        
+        if (ischar (varargin{pair_index + 1}))
+          if (! any (strcmp (lower (varargin{pair_index + 1}), ...
+                    {"sqeuclidean", "euclidean", "cityblock", "cosine", ...
+                     "correlation", "hamming", "jaccard"})))
+            error ("evalclusters: unknown distance criterion '%s'", ...
+                    varargin{pair_index + 1});
+          endif
+        elseif (! isa (varargin{pair_index + 1}, "function_handle") ||
+                ! ((isvector (varargin{pair_index + 1}) && ...
+                    isnumeric (varargin{pair_index + 1}))))
+          error ("evalclusters: invalid distance metric");
+        endif
+        distance = varargin{pair_index + 1};
+      case "clusterpriors"
+        ## used by silhouette evaluation
+        if (! strcmp (criterion, "silhouette"))
+          error (["evalclusters: cluster prior probabilities cannot be used "...
+                  "with '%s' criterion"], criterion);
+        endif
+        
+        if (any (strcmp (lower (varargin{pair_index + 1}), ...
+                                {"empirical", "equal"})))
+          clusterpriors = lower (varargin{pair_index + 1});
+        else
+          error ("evalclusters: invalid cluster prior probabilities value");
+        endif
+      case "b"
+        ## used by gap evaluation
+        if (! isnumeric (varargin{pair_index + 1}) || ...
+            ! isscalar (varargin{pair_index + 1}) || ...
+            varargin{pair_index + 1} != floor (varargin{pair_index + 1}) || ...
+            varargin{pair_index + 1} < 1)
+          error ("evalclusters: b must a be positive integer number");
+        endif
+        b = varargin{pair_index + 1};
+      case "referencedistribution"
+        ## used by gap evaluation
+        if (! ischar (varargin{pair_index + 1}) || any (strcmp (lower ...
+            (varargin{pair_index + 1}), {"pca", "uniform"})))
+          error (["evalclusters: the reference distribution must be either" ...
+                  "'PCA' or 'uniform'"]);
+        endif
+        referencedistribution = lower (varargin{pair_index + 1});
+      case "searchmethod"
+        ## used by gap evaluation
+        if (! ischar (varargin{pair_index + 1}) || any (strcmp (lower ...
+            (varargin{pair_index + 1}), {"globalmaxse", "uniform"})))
+          error (["evalclusters: the search method must be either" ...
+                  "'globalMaxSE' or 'firstmaxse'"]);
+        endif
+        searchmethod = lower (varargin{pair_index + 1});
+      otherwise
+        error ("evalclusters: unknown property %s", varargin{pair_index});
+    endswitch
+    
+    pair_index += 2;
+  endwhile 
+
+  ## check if there are parameters without a value or a name left
+  if (nargin - 2 - pair_index)
+    if (ischar (varargin{pair_index}))
+      error ("evalclusters: invalid parameter '%s'", varargin{pair_index});
+    else 
+      error ("evalclusters: invalid parameter '%d'", varargin{pair_index});
+    endif
+  endif 
+  
+  ## another check on klist
+  if (isempty (klist) && (ischar (clust) || isa (clust, "function_handle")))
+    error (["evalclusters: 'KList' can be empty only when 'clust' "...
+                    "is a matrix"]);
+  endif
+  
+  ## main
+  switch (lower (criterion))
+    case "calinskiharabasz"
+      ## further compatibility checks between the chosen parameters are 
+      ## delegated to the class constructor
+      if (isempty (klist))
+        klist = 1 : columns (clust);
+      endif
+      cc = CalinskiHarabaszEvaluation (x, clust, klist);
+    case "daviesbouldin"
+      ## further compatibility checks between the chosen parameters are 
+      ## delegated to the class constructor
+      if (isempty (klist))
+        klist = 1 : columns (clust);
+      endif
+      cc = DaviesBouldinEvaluation (x, clust, klist);
+    case "silhouette"
+      ## further compatibility checks between the chosen parameters are 
+      ## delegated to the class constructor
+      if (isempty (klist))
+        klist = 1 : columns (clust);
+      endif
+      cc = SilhouetteEvaluation (x, clust, klist, distance, clusterpriors);
+    case "gap"
+      ## gap cannot be used with a pre-computed solution, i.e. a matrix for
+      ## 'clust', and klist must be specified
+      if (isnumeric (clust))
+        error (["evalclusters: 'clust' must be a clustering algorithm when "... 
+                "using the gap criterion"]);
+      endif
+      if (isempty (klist))
+        error (["evalclusters: 'klist' cannot be empty when using the gap " ...
+                "criterion"]);
+      endif
+      cc = GapEvaluation (x, clust, klist, b, distance, ...
+                          referencedistribution, searchmethod);
+    otherwise
+      error ("evalclusters: invalid criterion '%s'", criterion);
+  endswitch
+endfunction
+
+## input tests
+%!error evalclusters ()
+%!error evalclusters ([1 1;0 1])
+%!error evalclusters ([1 1;0 1], "kmeans")
+%!error <'x' must be a numeric*> evalclusters ("abc", "kmeans", "gap")
+%!error <unknown clustering*> evalclusters ([1 1;0 1], "xxx", "gap")
+%!error <invalid matrix*> evalclusters ([1 1;0 1], [1 2], "gap")
+%!error <invalid argument*> evalclusters ([1 1;0 1], 1.2, "gap")
+%!error <invalid criterion*> evalclusters ([1 1;0 1], [1; 2], 123)
+%!error <unknown criterion*> evalclusters ([1 1;0 1], [1; 2], "xxx")
+%!error <'KList' can be empty*> evalclusters ([1 1;0 1], "kmeans", "gap")
+%!error <invalid parameter*> evalclusters ([1 1;0 1], [1; 2], "gap", 1)
+%!error <invalid property*> evalclusters ([1 1;0 1], [1; 2], "gap", 1, 1)
+%!error <unknown property*> evalclusters ([1 1;0 1], [1; 2], "gap", "xxx", 1)
+%!error <'KList'*> evalclusters ([1 1;0 1], [1; 2], "gap", "KList", [-1 0])
+%!error <'KList'*> evalclusters ([1 1;0 1], [1; 2], "gap", "KList", [1 .5])
+%!error <'KList'*> evalclusters ([1 1;0 1], [1; 2], "gap", "KList", [1 1; 1 1])
+%!error <unknown distance*> evalclusters ([1 1;0 1], [1; 2], "gap", ...
+%!                                        "distance", "a")
+%!error <distance metric*> evalclusters ([1 1;0 1], [1; 2], "daviesbouldin", ...
+%!                                       "distance", "a")
+%!error <cluster prior*> evalclusters ([1 1;0 1], [1; 2], "gap", ...
+%!                                     "clusterpriors", "equal")
+%!error <invalid cluster prior*> evalclusters ([1 1;0 1], [1; 2], ...
+%!                                         "silhouette", "clusterpriors", "xxx")
+%!error <'clust' must be a clustering*> evalclusters ([1 1;0 1], [1; 2], "gap")
+
+## demonstration
+%!demo
+%! load fisheriris;
+%! eva = evalclusters(meas, "kmeans", "calinskiharabasz", "KList", [1:6])
+%! figure();
+%! eva.plot;
\ No newline at end of file
diff -r d354b1b49cad -r 933bc06e11e4 inst/private/CalinskiHarabaszEvaluation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/private/CalinskiHarabaszEvaluation.m	Sat Apr 03 17:51:12 2021 +0200
@@ -0,0 +1,202 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+classdef CalinskiHarabaszEvaluation < ClusterCriterion
+  ## -*- texinfo -*-
+  ## @deftypefn {Function File} {@var{eva} =} evalclusters (@var{x}, @var{clust}, @qcode{CalinskiHarabasz})
+  ## @deftypefnx {Function File} {@var{eva} =} evalclusters (@dots{}, @qcode{Name}, @qcode{Value})
+  ##
+  ## A Calinski-Harabasz object to evaluate clustering solutions.
+  ##
+  ## A @code{CalinskiHarabaszEvaluation} object is a @code{ClusterCriterion}
+  ## object used to evaluate clustering solutions using the Calinski-Harabasz
+  ## criterion.
+  ##
+  ## The Calinski-Harabasz index is based on the ratio between SSb and SSw.
+  ## SSb is the overall variance between clusters, that is the variance of the
+  ## distances between the centroids.
+  ## SSw is the overall variance within clusters, that is the sum of the
+  ## variances of the distances between each datapoint and its centroid.
+  ##
+  ## The best solution according to the Calinski-Harabasz criterion is the one
+  ## that scores the highest value.
+  ## @end deftypefn
+  ##
+  ## @seealso{ClusterCriterion, evalclusters}
+
+  properties (GetAccess = public, SetAccess = private)
+
+  endproperties
+  
+  properties (Access = protected)
+    Centroids = {}; # a list of the centroids for every solution
+  endproperties
+  
+  methods (Access = public)
+    ## constructor
+    function this = CalinskiHarabaszEvaluation (x, clust, KList)
+      this@ClusterCriterion(x, clust, KList);
+      
+      this.CriterionName = "CalinskiHarabasz";
+      this.evaluate(this.InspectedK); # evaluate the list of cluster numbers
+    endfunction
+    
+    ## set functions
+
+    ## addK
+    ## add new cluster sizes to evaluate
+    function this = addK (this, K)
+      addK@ClusterCriterion(this, K);
+      
+      ## if we have new data, we need a new evaluation
+      if (this.OptimalK == 0)
+        Centroids_tmp = {};
+        pS = 0; # position shift of the elements of Centroids 
+        for iter = 1 : length (this.InspectedK)
+          ## reorganize Centroids according to the new list of cluster numbers
+          if (any (this.InspectedK(iter) == K))
+            pS += 1;
+          else
+            Centroids_tmp{iter} = this.Centroids{iter - pS};
+          endif
+        endfor
+        this.Centroids = Centroids_tmp;
+        this.evaluate(K); # evaluate just the new cluster numbers
+      endif
+    endfunction
+    
+    ## compact
+    ## ...
+    function this = compact (this)
+      # FIXME: stub!
+      warning ("CalinskiHarabaszEvaluation: compact is unavailable");
+    endfunction
+  endmethods
+  
+  methods (Access = protected)
+    ## evaluate
+    ## do the evaluation
+    function this = evaluate (this, K)
+      ## use complete observations only
+      UsableX = this.X(find (this.Missing == false), :);
+      if (! isempty (this.ClusteringFunction))
+        ## build the clusters
+        for iter = 1 : length (this.InspectedK)
+          ## do it only for the specified K values
+          if (any (this.InspectedK(iter) == K))
+            if (isa (this.ClusteringFunction, "function_handle"))
+              ## custom function
+              ClusteringSolution = ...
+                this.ClusteringFunction(UsableX, this.InspectedK(iter));
+              if (ismatrix (ClusteringSolution) && ...
+                  rows (ClusteringSolution) == this.NumObservations && ...
+                  columns (ClusteringSolution) == this.P)
+                ## the custom function returned a matrix:
+                ## we take the index of the maximum value for every row
+                [~, this.ClusteringSolutions(:, iter)] = ...
+                  max (ClusteringSolution, [], 2);
+              elseif (iscolumn (ClusteringSolution) &&
+                      length (ClusteringSolution) == this.NumObservations)
+                this.ClusteringSolutions(:, iter) = ClusteringSolution;
+              elseif (isrow (ClusteringSolution) &&
+                      length (ClusteringSolution) == this.NumObservations)
+                this.ClusteringSolutions(:, iter) = ClusteringSolution';
+              else
+                error (["CalinskiHarabaszEvaluation: invalid return value "...
+                        "from custom clustering function"]);
+              endif
+              this.ClusteringSolutions(:, iter) = ...
+                this.ClusteringFunction(UsableX, this.InspectedK(iter));
+            else
+              switch (this.ClusteringFunction)
+                case "kmeans"
+                  [this.ClusteringSolutions(:, iter), this.Centroids{iter}] =...
+                    kmeans (UsableX, this.InspectedK(iter),  ...
+                    "Distance", "sqeuclidean", "EmptyAction", "singleton", ...
+                    "Replicates", 5);
+                case "linkage"
+                  ## use clusterdata
+                  this.ClusteringSolutions(:, iter) = clusterdata (UsableX, ...
+                    "MaxClust", this.InspectedK(iter), ...
+                    "Distance", "euclidean", "Linkage", "ward");
+                  this.Centroids{iter} = this.computeCentroids (UsableX, iter);
+                case "gmdistribution"
+                  gmm = fitgmdist (UsableX, this.InspectedK(iter), ...
+                        "SharedCov", true, "Replicates", 5);
+                  this.ClusteringSolutions(:, iter) = cluster (gmm, UsableX);
+                  this.Centroids{iter} = gmm.mu;
+                otherwise
+                  error (["CalinskiHarabaszEvaluation: unexpected error, " ...
+                          "report this bug"]);
+              endswitch
+            endif
+          endif
+        endfor
+      endif
+      
+      ## get the criterion values for every clustering solution
+      for iter = 1 : length (this.InspectedK)
+        ## do it only for the specified K values
+        if (any (this.InspectedK(iter) == K))
+          ## not defined for one cluster
+          if (this.InspectedK(iter) == 1)
+            this.CriterionValues(iter) = NaN;
+            continue;
+          endif
+          
+          ## Caliński-Harabasz index
+          ## reference: calinhara function from the fpc package of R,
+          ##            by Christian Hennig
+          ##            https://CRAN.R-project.org/package=fpc
+          W = zeros (columns (UsableX)); # between clusters covariance
+          for i = 1 : this.InspectedK(iter)
+            vIndicesI = find (this.ClusteringSolutions(:, iter) == i);
+            ni = length (vIndicesI); # size of cluster i
+            if (ni == 1)
+              ## if the cluster has just one member the covariance is zero
+              continue;
+            endif
+            ## weighted update of the covariance matrix
+            W += cov (UsableX(vIndicesI, :)) * (ni - 1);
+          endfor
+          S = (this.NumObservations - 1) * cov (UsableX); # within clusters cov.
+          B = S - W; # between clusters means
+          
+          ## tr(B) / tr(W) * (N-k) / (k-1)
+          this.CriterionValues(iter) = (this.NumObservations - ...
+            this.InspectedK(iter)) * trace (B) / ...
+            ((this.InspectedK(iter) - 1) * trace (W));
+        endif
+      endfor
+      
+      [~, this.OptimalIndex] = max (this.CriterionValues);
+      this.OptimalK = this.InspectedK(this.OptimalIndex(1));
+      this.OptimalY = this.ClusteringSolutions(:, this.OptimalIndex(1));
+    endfunction
+  endmethods
+  
+  methods (Access = private)
+    ## computeCentroids
+    ## compute the centroids if they are not available by other means
+    function C = computeCentroids (this, X, index)
+      C = zeros (this.InspectedK(index), columns (X));
+      
+      for iter = 1 : this.InspectedK(index)
+        vIndicesI = find (this.ClusteringSolutions(:, index) == iter);
+        C(iter, :) = mean (X(vIndicesI, :));
+      endfor
+    endfunction
+  endmethods
+endclassdef
diff -r d354b1b49cad -r 933bc06e11e4 inst/private/ClusterCriterion.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/private/ClusterCriterion.m	Sat Apr 03 17:51:12 2021 +0200
@@ -0,0 +1,240 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+classdef ClusterCriterion < handle
+  ## -*- texinfo -*-
+  ## @deftypefn {} {} ClusterCriterion ()
+  ##
+  ## A clustering evaluation object as created by @code{evalclusters}.
+  ##
+  ## @code{ClusterCriterion} is a superclass for clustering evaluation objects
+  ## as created by @code{evalclusters}.
+  ##
+  ## List of public properties:
+  ## @table @code
+  ## @item @qcode{ClusteringFunction} 
+  ## a valid clustering funtion name or function handle. It can be empty if the
+  ## clustering solutions are passed as an input matric.
+  ##
+  ## @item @qcode{CriterionName} 
+  ## a valid criterion name to evaluate the clustering solutions.
+  ##
+  ## @item @qcode{CriterionValues} 
+  ## a vector of values as generated by the evaluation criterion for each
+  ## clustering solution.
+  ##
+  ## @item @qcode{InspectedK} 
+  ## the list of proposed cluster numbers.
+  ##
+  ## @item @qcode{Missing} 
+  ## a logical vector of missing observations. When there are @code{NaN} values
+  ## in the data matrix, the corresponding observation is excluded.
+  ##
+  ## @item @qcode{NumObservations} 
+  ## the number of non-missing observations in the data matrix.
+  ##
+  ## @item @qcode{OptimalK} 
+  ## the optimal number of clusters.
+  ##
+  ## @item @qcode{OptimalY} 
+  ## the clustering solution corresponding to @code{OptimalK} .
+  ##
+  ## @item @qcode{X} 
+  ## the data matrix.
+  ##
+  ## @end table
+  ##
+  ## List of public methods:
+  ## @table @code
+  ## @item @qcode{addK} 
+  ## add a list of numbers of clusters to evaluate.
+  ##
+  ## @item @qcode{compact} 
+  ## return a compact clustering evaluation object. Not implemented
+  ##
+  ## @item @qcode{plot} 
+  ## plot the clustering evaluation values against the corresponding number of
+  ## clusters.
+  ##
+  ## @end table
+  ## @end deftypefn
+  ##
+  ## @seealso{CalinskiHarabaszEvaluation, DaviesBouldinEvaluation, evalclusters,
+  ## GapEvaluation, SilhouetteEvaluation}
+  
+  properties (Access = public)
+    ## public properties
+  endproperties
+
+  properties (GetAccess = public, SetAccess = protected)
+    ClusteringFunction = "";
+    CriterionName = "";
+    CriterionValues = [];
+    InspectedK = [];
+    Missing = [];
+    NumObservations = 0;
+    OptimalK = 0;
+    OptimalY = [];
+    X = [];
+  endproperties
+
+  properties (Access = protected)
+    N = 0; # number of observations
+    P = 0; # number of variables
+    ClusteringSolutions = []; #
+    OptimalIndex = 0; # index of the optimal K
+  endproperties
+  
+  methods (Access = public)
+    ## constructor
+    function this = ClusterCriterion (x, clust, KList)
+      ## parsing input data
+      if ((! ismatrix (x)) || (! isnumeric (x)))
+        error ("ClusterCriterion: 'x' must be a numeric matrix");
+      endif
+      this.X = x;
+      this.N = rows (this.X);
+      this.P = columns (this.X);
+      ## look for missing values
+      for iter = 1 : this.N
+        if (any (find (x(iter, :) == NaN)))
+          this.Missing(iter) = true;
+        else
+          this.Missing(iter) = false;
+        endif
+      endfor
+      ## number of usable observations
+      this.NumObservations = sum (this.Missing == false);
+      
+      ## parsing the clustering algorithm
+      if (ischar (clust))
+        if (any (strcmp (lower (clust), {"kmeans", "linkage", ...
+                                         "gmdistribution"})))
+          this.ClusteringFunction = lower (clust);
+        else
+          error ("ClusterCriterion: unknown clustering algorithm '%s'", clust);
+        endif
+      elseif (isa (clust, "function_handle"))
+        this.ClusteringFunction = clust;
+      elseif (ismatrix (clust))
+        if (isnumeric (clust)  && (length (size (clust)) == 2) && ...
+            (rows (clust) == this.N))
+          this.ClusteringFunction = "";
+          this.ClusteringSolutions = clust(find (this.Missing == false), :);
+        else
+          error ("ClusterCriterion: invalid matrix of clustering solutions");
+        endif
+      else
+        error ("ClusterCriterion: invalid argument");
+      endif
+      
+      ## parsing the list of cluster sizes to inspect
+      this.InspectedK = parseKList (this, KList);
+    endfunction
+    
+    ## addK
+    ## add one or more new cluster sizes to evaluate
+    function this = addK (this, k)
+      ## -*- texinfo -*-
+      ## @deftypefn {} {} ClusterCriterion.addK (@var{K})
+      ##
+      ## Add an array of cluster numbers to inspect to the evaluation object.
+      ## @end deftypefn
+      
+      ## if there is not a clustering function, then we are using a predefined
+      ## set of clustering solutions, hence we cannot redefine the number of
+      ## solutions
+      if (isempty (this.ClusteringFunction))
+        warning (["ClusterCriterion: cannot redefine the list of cluster"...
+          "numbers to evaluate when there is not a clustering function"]);
+        return;
+      endif
+      
+      ## otherwise go on
+      newList = this.parseKList ([this.InspectedK k]);
+      
+      ## check if the list has changed
+      if (length (newList) == length (this.InspectedK))
+        warning ("ClusterCriterion: the list has not changed");
+      else
+        ## update ClusteringSolutions and CriterionValues
+        ClusteringSolutions_tmp = zeros (this.NumObservations, ...
+                                    length (newList));
+        CriterionValues_tmp = zeros (length (newList), 1);
+        for iter = 1 : length (this.InspectedK)
+          idx = find (newList == this.InspectedK(iter));
+          
+          if (! isempty (idx))
+            ClusteringSolutions_tmp(:, idx) = this.ClusteringSolutions(:, iter);
+            CriterionValues_tmp(idx) = this.CriterionValues(iter);
+          endif
+        endfor
+        this.ClusteringSolutions = ClusteringSolutions_tmp;
+        this.CriterionValues = CriterionValues_tmp;
+        
+        ## reset the old results
+        this.OptimalK = 0;
+        this.OptimalY = [];
+        this.OptimalIndex = 0;
+        
+        ## update the list of cluster numbers to evaluate
+        this.InspectedK = newList;
+      endif
+    endfunction
+    
+    ## plot
+    ## plot the CriterionValues against InspectedK and return a handle to the
+    ## plot
+    function h = plot (this)
+      ## -*- texinfo -*-
+      ## @deftypefn {} {} ClusterCriterion.plot ()
+      ##
+      ## Plot the evaluation results.
+      ## @end deftypefn
+      
+      yLabel = sprintf ("%s value", this.CriterionName);
+      h = gca ();
+      hold on;
+      plot (this.InspectedK, this.CriterionValues, "bo-");
+      plot (this.OptimalK, this.CriterionValues(this.OptimalIndex), "b*");
+      xlabel ("number of clusters");
+      ylabel (yLabel);
+      hold off;
+    endfunction
+  endmethods
+  
+  methods (Abstract = true)
+    function compact ()
+      ## -*- texinfo -*-
+      ## @deftypefn {} {@var{eva} =} ClusterCriterion.compact ()
+      ##
+      ## Return a compact evaluation object.
+      ## @end deftypefn
+    endfunction
+  endmethods
+  
+  methods (Access = private)
+    ## check if a list of cluster sizes is correct
+    function retList = parseKList (this, KList)
+      if (isnumeric (KList) && isvector (KList) && all (find (KList > 0)) && ...
+          all (floor (KList) == KList))
+        retList = unique (KList);
+      else
+        error (["ClusterCriterion: the list of cluster sizes must be an " ...
+                "array of positive integer numbers"]);
+      endif
+    endfunction
+  endmethods
+endclassdef
diff -r d354b1b49cad -r 933bc06e11e4 inst/private/DaviesBouldinEvaluation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/private/DaviesBouldinEvaluation.m	Sat Apr 03 17:51:12 2021 +0200
@@ -0,0 +1,205 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+classdef DaviesBouldinEvaluation < ClusterCriterion
+  ## -*- texinfo -*-
+  ## @deftypefn {Function File} {@var{eva} =} evalclusters (@var{x}, @var{clust}, @qcode{DaviesBouldin})
+  ## @deftypefnx {Function File} {@var{eva} =} evalclusters (@dots{}, @qcode{Name}, @qcode{Value})
+  ##
+  ## A Davies-Bouldin object to evaluate clustering solutions.
+  ##
+  ## A @code{DaviesBouldinEvaluation} object is a @code{ClusterCriterion}
+  ## object used to evaluate clustering solutions using the Davies-Bouldin
+  ## criterion.
+  ##
+  ## The Davies-Bouldin criterion is based on the ratio between the distances
+  ## between clusters and within clusters, that is the distances between the 
+  ## centroids and the distances between each datapoint and its centroid.
+  ##
+  ## The best solution according to the Davies-Bouldin criterion is the one
+  ## that scores the lowest value.
+  ## @end deftypefn
+  ##
+  ## @seealso{ClusterCriterion, evalclusters}
+  
+  properties (GetAccess = public, SetAccess = private)
+
+  endproperties
+  
+  properties (Access = protected)
+    Centroids = {}; # a list of the centroids for every solution
+  endproperties
+  
+  methods (Access = public)
+    ## constructor
+    function this = DaviesBouldinEvaluation (x, clust, KList)
+      this@ClusterCriterion(x, clust, KList);
+      
+      this.CriterionName = "DaviesBouldin";
+      this.evaluate(this.InspectedK); # evaluate the list of cluster numbers
+    endfunction
+    
+    ## set functions
+
+    ## addK
+    ## add new cluster sizes to evaluate
+    function this = addK (this, K)
+      addK@ClusterCriterion(this, K);
+      
+      ## if we have new data, we need a new evaluation
+      if (this.OptimalK == 0)
+        Centroids_tmp = {};
+        pS = 0; # position shift of the elements of Centroids 
+        for iter = 1 : length (this.InspectedK)
+          ## reorganize Centroids according to the new list of cluster numbers
+          if (any (this.InspectedK(iter) == K))
+            pS += 1;
+          else
+            Centroids_tmp{iter} = this.Centroids{iter - pS};
+          endif
+        endfor
+        this.Centroids = Centroids_tmp;
+        this.evaluate(K); # evaluate just the new cluster numbers
+      endif
+    endfunction
+    
+    ## compact
+    ## ...
+    function this = compact (this)
+      # FIXME: stub!
+      warning ("DaviesBouldinEvaluation: compact is unavailable");
+    endfunction
+  endmethods
+  
+  methods (Access = protected)
+    ## evaluate
+    ## do the evaluation
+    function this = evaluate (this, K)
+      ## use complete observations only
+      UsableX = this.X(find (this.Missing == false), :);
+      if (! isempty (this.ClusteringFunction))
+        ## build the clusters
+        for iter = 1 : length (this.InspectedK)
+          ## do it only for the specified K values
+          if (any (this.InspectedK(iter) == K))
+            if (isa (this.ClusteringFunction, "function_handle"))
+              ## custom function
+              ClusteringSolution = ...
+                this.ClusteringFunction(UsableX, this.InspectedK(iter));
+              if (ismatrix (ClusteringSolution) && ...
+                  rows (ClusteringSolution) == this.NumObservations && ...
+                  columns (ClusteringSolution) == this.P)
+                ## the custom function returned a matrix:
+                ## we take the index of the maximum value for every row
+                [~, this.ClusteringSolutions(:, iter)] = ...
+                  max (ClusteringSolution, [], 2);
+              elseif (iscolumn (ClusteringSolution) &&
+                      length (ClusteringSolution) == this.NumObservations)
+                this.ClusteringSolutions(:, iter) = ClusteringSolution;
+              elseif (isrow (ClusteringSolution) &&
+                      length (ClusteringSolution) == this.NumObservations)
+                this.ClusteringSolutions(:, iter) = ClusteringSolution';
+              else
+                error (["DaviesBouldinEvaluation: invalid return value "...
+                        "from custom clustering function"]);
+              endif
+              this.ClusteringSolutions(:, iter) = ...
+                this.ClusteringFunction(UsableX, this.InspectedK(iter));
+            else
+              switch (this.ClusteringFunction)
+                case "kmeans"
+                  [this.ClusteringSolutions(:, iter), this.Centroids{iter}] =...
+                    kmeans (UsableX, this.InspectedK(iter),  ...
+                    "Distance", "sqeuclidean", "EmptyAction", "singleton", ...
+                    "Replicates", 5);
+                case "linkage"
+                  ## use clusterdata
+                  this.ClusteringSolutions(:, iter) = clusterdata (UsableX, ...
+                    "MaxClust", this.InspectedK(iter), ...
+                    "Distance", "euclidean", "Linkage", "ward");
+                  this.Centroids{iter} = this.computeCentroids (UsableX, iter);
+                case "gmdistribution"
+                  gmm = fitgmdist (UsableX, this.InspectedK(iter), ...
+                        "SharedCov", true, "Replicates", 5);
+                  this.ClusteringSolutions(:, iter) = cluster (gmm, UsableX);
+                  this.Centroids{iter} = gmm.mu;
+                otherwise
+                  error (["DaviesBouldinEvaluation: unexpected error, " ...
+                          "report this bug"]);
+              endswitch
+            endif
+          endif
+        endfor
+      endif
+      
+      ## get the criterion values for every clustering solution
+      for iter = 1 : length (this.InspectedK)
+        ## do it only for the specified K values
+        if (any (this.InspectedK(iter) == K))
+          ## not defined for one cluster
+          if (this.InspectedK(iter) == 1)
+            this.CriterionValues(iter) = NaN;
+            continue;
+          endif
+          
+          ## Davies-Bouldin value
+          ## an evaluation of the ratio between within-cluster and 
+          ## between-cluster distances
+          
+          ## mean distances between cluster members and their centroid
+          vD = zeros (this.InspectedK(iter), 1);
+          for i = 1 : this.InspectedK(iter)
+            vIndicesI = find (this.ClusteringSolutions(:, iter) == i);
+            vD(i) = mean (vecnorm (UsableX(vIndicesI, :) - ...
+                                   this.Centroids{iter}(i, :), 2, 2));
+          endfor
+          
+          ## within-to-between cluster distance ratio
+          Dij = zeros (this.InspectedK(iter));
+          for i = 1 : (this.InspectedK(iter) - 1)
+            for j = (i + 1) : this.InspectedK(iter)
+              ## centroid to centroid distance
+              dij = vecnorm (this.Centroids{iter}(i, :) - ...
+                             this.Centroids{iter}(j, :));
+              ## within-to-between cluster distance ratio for clusters i and j
+              Dij(i, j) = (vD(i) + vD(j)) / dij;
+            endfor
+          endfor
+          
+          ## ( max_j D1j + max_j D2j + ... + max_j Dkj) / k
+          this.CriterionValues(iter) = sum (max (Dij(i, :), [], 2)) / ...
+                                            this.InspectedK(iter);
+        endif
+      endfor
+      
+      [~, this.OptimalIndex] = min (this.CriterionValues);
+      this.OptimalK = this.InspectedK(this.OptimalIndex(1));
+      this.OptimalY = this.ClusteringSolutions(:, this.OptimalIndex(1));
+    endfunction
+  endmethods
+  
+  methods (Access = private)
+    ## computeCentroids
+    ## compute the centroids if they are not available by other means
+    function C = computeCentroids (this, X, index)
+      C = zeros (this.InspectedK(index), columns (X));
+      
+      for iter = 1 : this.InspectedK(index)
+        vIndicesI = find (this.ClusteringSolutions(:, index) == iter);
+        C(iter, :) = mean (X(vIndicesI, :));
+      endfor
+    endfunction
+  endmethods
+endclassdef
diff -r d354b1b49cad -r 933bc06e11e4 inst/private/GapEvaluation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/private/GapEvaluation.m	Sat Apr 03 17:51:12 2021 +0200
@@ -0,0 +1,362 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+classdef GapEvaluation < ClusterCriterion
+  ## -*- texinfo -*-
+  ## @deftypefn {Function File} {@var{eva} =} evalclusters (@var{x}, @var{clust}, @qcode{gap})
+  ## @deftypefnx {Function File} {@var{eva} =} evalclusters (@dots{}, @qcode{Name}, @qcode{Value})
+  ##
+  ## A gap object to evaluate clustering solutions.
+  ##
+  ## A @code{GapEvaluation} object is a @code{ClusterCriterion}
+  ## object used to evaluate clustering solutions using the gap criterion,
+  ## which is a mathematical formalization of the elbow method.
+  ##
+  ## List of public properties specific to @code{SilhouetteEvaluation}:
+  ## @table @code
+  ## @item @qcode{B} 
+  ## the number of reference datasets to generate.
+  ##
+  ## @item @qcode{Distance} 
+  ## a valid distance metric name, or a function handle as accepted by the
+  ## @code{pdist} function.
+  ##  
+  ## @item @qcode{ExpectedLogW} 
+  ## a vector of the expected values for the logarithm of the within clusters 
+  ## dispersion.
+  ##  
+  ## @item @qcode{LogW} 
+  ## a vector of the values of the logarithm of the within clusters dispersion.
+  ##
+  ## @item @qcode{ReferenceDistribution} 
+  ## a valid name for the reference distribution, namely: @code{PCA} (default)
+  ## or @code{uniform}.
+  ##  
+  ## @item @qcode{SE} 
+  ## a vector of the standard error of the expected values for the logarithm 
+  ## of the within clusters dispersion.
+  ##  
+  ## @item @qcode{SearchMethod} 
+  ## a valid name for the search method to use: @code{globalMaxSE} (default) or
+  ## @code{firstMaxSE}.
+  ##  
+  ## @item @qcode{StdLogW} 
+  ## a vector of the standard deviation of the expected values for the logarithm 
+  ## of the within clusters dispersion.
+  ## @end table
+  ## 
+  ## The best solution according to the gap criterion depends on the chosen
+  ## search method. When the search method is @code{globalMaxSE}, the chosen
+  ## gap value is the smaller one which is inside a standard error from the
+  ## max gap value; when the search method is @code{firstMaxSE}, the chosen
+  ## gap value is the first one which is inside a standard error from the next
+  ## gap value.
+  ## @end deftypefn
+  ##
+  ## @seealso{ClusterCriterion, evalclusters}
+  
+  properties (GetAccess = public, SetAccess = private)
+    B = 0; # number of reference datasets
+    Distance = ""; # pdist parameter
+    ReferenceDistribution = ""; # distribution to use as reference
+    SearchMethod = ""; # the method do identify the optimal number of clusters
+    ExpectedLogW = []; # expected value for the natural logarithm of W
+    LogW = []; # natural logarithm of W
+    SE = []; # standard error for the natural logarithm of W
+    StdLogW = []; # standard deviation of the natural logarithm of W
+  endproperties
+  
+  properties (Access = protected)
+    DistanceVector = []; # vector of pdist distances
+    mExpectedLogW = []; # the result of the Monte-Carlo simulations
+  endproperties
+  
+  methods (Access = public)
+    ## constructor
+    function this = GapEvaluation (x, clust, KList, b = 100, ...
+                    distanceMetric = "sqeuclidean", ...
+                    referenceDistribution = "pca", searchMethod = "globalmaxse")
+      this@ClusterCriterion(x, clust, KList);
+      
+      ## parsing the distance criterion
+      if (ischar (distanceMetric))
+        if (any (strcmp (lower (distanceMetric), {"sqeuclidean", ...
+                  "euclidean", "cityblock", "cosine", "correlation", ...
+                  "hamming", "jaccard"})))
+          this.Distance = lower (distanceMetric);
+          
+          ## kmeans can use only a subset
+          if (strcmp (clust, "kmeans") && any (strcmp (this.Distance, ...
+              {"euclidean", "jaccard"})))
+            error (["GapEvaluation: invalid distance criterion '%s' "...
+                    "for 'kmeans'"], distanceMetric);
+          endif
+        else
+          error ("GapEvaluation: unknown distance criterion '%s'", ...
+                 distanceMetric);
+        endif
+      elseif (isa (distanceMetric, "function_handle"))
+        this.Distance = distanceMetric;
+        
+        ## kmeans cannot use a function handle
+        if (strcmp (clust, "kmeans"))
+          error ("GapEvaluation: invalid distance criterion for 'kmeans'");
+        endif
+      elseif (isvector (distanceMetric) && isnumeric (distanceMetric))
+        this.Distance = "";
+        this.DistanceVector = distanceMetric; # the validity check is delegated
+        
+        ## kmeans cannot use a distance vector
+        if (strcmp (clust, "kmeans"))
+          error (["GapEvaluation: invalid distance criterion for "...
+                  "'kmeans'"]);
+        endif
+      else
+        error ("GapEvaluation: invalid distance metric");
+      endif
+
+      ## B: number of Monte-Carlo iterations
+      if (! isnumeric (b) || ! isscalar (b) || b != floor (b) || b < 1)
+        error ("GapEvaluation: b must a be positive integer number");
+      endif
+      this.B = b;
+      
+      ## reference distribution
+      if (! ischar (referenceDistribution) || ! any (strcmp (lower ...
+          (referenceDistribution), {"pca", "uniform"})))
+        error (["GapEvaluation: the reference distribution must be either" ...
+                "'PCA' or 'uniform'"]);
+      elseif (strcmp (lower (referenceDistribution), "pca"))
+        warning (["GapEvaluation: 'PCA' distribution not implemented, " ...
+                  "using 'uniform'"]);
+      endif
+      this.ReferenceDistribution = lower (referenceDistribution);
+      
+      if (! ischar (searchMethod) || ! any (strcmp (lower ...
+          (searchMethod), {"globalmaxse", "firstmaxse"})))
+        error (["evalclusters: the search method must be either" ...
+                "'globalMaxSE' or 'firstMaxSE'"]);
+      endif
+      this.SearchMethod = lower (searchMethod);
+      
+      ## a matrix to store the results from the Monte-Carlo runs
+      this.mExpectedLogW = zeros (this.B, length (this.InspectedK));
+      
+      this.CriterionName = "gap";
+      this.evaluate(this.InspectedK); # evaluate the list of cluster numbers
+    endfunction
+    
+    ## set functions
+
+    ## addK
+    ## add new cluster sizes to evaluate
+    function this = addK (this, K)
+      addK@ClusterCriterion(this, K);
+      
+      ## if we have new data, we need a new evaluation
+      if (this.OptimalK == 0)
+        mExpectedLogW_tmp = zeros (this.B, length (this.InspectedK));
+        pS = 0; # position shift
+        for iter = 1 : length (this.InspectedK)
+          ## reorganize all the arrays according to the new list
+          ## of cluster numbers
+          if (any (this.InspectedK(iter) == K))
+            pS += 1;
+          else
+            mExpectedLogW_tmp(:, iter) = this.mExpectedLogW(:, iter - pS);
+          endif
+        endfor
+        this.mExpectedLogW = mExpectedLogW_tmp;
+        
+        this.evaluate(K); # evaluate just the new cluster numbers
+      endif
+    endfunction
+    
+    ## compact
+    ## ...
+    function this = compact (this)
+      # FIXME: stub!
+      warning ("GapEvaluation: compact is unavailable");
+    endfunction
+    
+    ## plot
+    ## plot the CriterionValues against InspectedK, show the standard deviation
+    ## and return a handle to the plot
+    function h = plot (this)
+      yLabel = sprintf ("%s value", this.CriterionName);
+      h = gca ();
+      hold on;
+      errorbar (this.InspectedK, this.CriterionValues, this.StdLogW);
+      plot (this.InspectedK, this.CriterionValues, "bo");
+      plot (this.OptimalK, this.CriterionValues(this.OptimalIndex), "b*");
+      xlabel ("number of clusters");
+      ylabel (yLabel);
+      hold off;
+    endfunction
+  endmethods
+  
+  methods (Access = protected)
+    ## evaluate
+    ## do the evaluation
+    function this = evaluate (this, K)      
+      ## Monte-Carlo runs
+      for mcrun = 1 : (this.B + 1)
+        ## use complete observations only
+        UsableX = this.X(find (this.Missing == false), :);
+        
+        ## the last run use tha actual data,
+        ## the others are Monte-Carlo runs with reconstructed data
+        if (mcrun <= this.B)
+          ## uniform distribution
+          colMins = min (UsableX);
+          colMaxs = max (UsableX);
+          for col = 1 : columns (UsableX)
+            UsableX(:, col) = colMins(col) + rand (this.NumObservations, 1) *...
+                              (colMaxs(col) - colMins(col));
+          endfor
+        endif  
+      
+        if (! isempty (this.ClusteringFunction))
+          ## build the clusters
+          for iter = 1 : length (this.InspectedK)
+            ## do it only for the specified K values
+            if (any (this.InspectedK(iter) == K))
+              if (isa (this.ClusteringFunction, "function_handle"))
+                ## custom function
+                ClusteringSolution = ...
+                  this.ClusteringFunction(UsableX, this.InspectedK(iter));
+                if (ismatrix (ClusteringSolution) && ...
+                    rows (ClusteringSolution) == this.NumObservations && ...
+                    columns (ClusteringSolution) == this.P)
+                  ## the custom function returned a matrix:
+                  ## we take the index of the maximum value for every row
+                  [~, this.ClusteringSolutions(:, iter)] = ...
+                    max (ClusteringSolution, [], 2);
+                elseif (iscolumn (ClusteringSolution) &&
+                        length (ClusteringSolution) == this.NumObservations)
+                  this.ClusteringSolutions(:, iter) = ClusteringSolution;
+                elseif (isrow (ClusteringSolution) &&
+                        length (ClusteringSolution) == this.NumObservations)
+                  this.ClusteringSolutions(:, iter) = ClusteringSolution';
+                else
+                  error (["GapEvaluation: invalid return value from " ...
+                          "custom clustering function"]);
+                endif
+                this.ClusteringSolutions(:, iter) = ...
+                  this.ClusteringFunction(UsableX, this.InspectedK(iter));
+              else
+                switch (this.ClusteringFunction)
+                  case "kmeans"
+                    this.ClusteringSolutions(:, iter) = kmeans (UsableX, ...
+                      this.InspectedK(iter), "Distance", this.Distance, ...
+                      "EmptyAction", "singleton", "Replicates", 5);
+                  case "linkage"
+                    if (! isempty (this.Distance))
+                      ## use clusterdata
+                      Distance_tmp = this.Distance;
+                      LinkageMethod = "average"; # for non euclidean methods
+                      if (strcmp (this.Distance, "sqeuclidean"))
+                        ## pdist uses different names for its algorithms
+                        Distance_tmp = "squaredeuclidean";
+                        LinkageMethod = "ward";
+                      elseif (strcmp (this.Distance, "euclidean"))
+                        LinkageMethod = "ward";
+                      endif
+                      this.ClusteringSolutions(:, iter) = clusterdata ...
+                        (UsableX, "MaxClust", this.InspectedK(iter), ...
+                        "Distance", Distance_tmp, "Linkage", LinkageMethod);
+                    else
+                      ## use linkage
+                      Z = linkage (this.DistanceVector, "average");
+                      this.ClusteringSolutions(:, iter) = ...
+                           cluster (Z, "MaxClust", this.InspectedK(iter));
+                    endif
+                  case "gmdistribution"
+                    gmm = fitgmdist (UsableX, this.InspectedK(iter), ...
+                          "SharedCov", true, "Replicates", 5);
+                    this.ClusteringSolutions(:, iter) = cluster (gmm, UsableX);
+                  otherwise
+                    ## this should not happen
+                    error (["GapEvaluation: unexpected error, " ...
+                           "report this bug"]);
+                endswitch
+              endif
+            endif
+          endfor
+        endif
+        
+        ## get the gap values for every clustering
+        distance_pdist = this.Distance;
+        if (strcmp (distance_pdist, "sqeuclidean"))
+          distance_pdist = "squaredeuclidean";
+        endif
+        
+        ## compute LogW
+        for iter = 1 : length (this.InspectedK)
+          ## do it only for the specified K values
+          if (any (this.InspectedK(iter) == K))
+            wk = 0;
+            for r = 1 : this.InspectedK(iter)
+              vIndicesR = find (this.ClusteringSolutions(:, iter) == r);
+              nr = length (vIndicesR);
+              Dr = pdist (UsableX(vIndicesR, :), distance_pdist);
+              wk += sum (Dr) / (2 * nr);
+            endfor
+            if (mcrun <= this.B)
+              this.mExpectedLogW(mcrun, iter) = log (wk);
+            else
+              this.LogW(iter) = log (wk);
+            endif
+          endif
+        endfor
+      endfor
+      
+      this.ExpectedLogW = mean (this.mExpectedLogW);
+      this.SE = sqrt ((1 + 1 / this.B) * sumsq (this.mExpectedLogW - ...
+                                               this.ExpectedLogW) / this.B);
+      this.StdLogW = std (this.mExpectedLogW);
+      this.CriterionValues = this.ExpectedLogW - this.LogW;
+      
+      this.OptimalIndex = this.gapSearch ();
+      this.OptimalK = this.InspectedK(this.OptimalIndex(1));
+      this.OptimalY = this.ClusteringSolutions(:, this.OptimalIndex(1));
+    endfunction
+    
+    ## gapSearch
+    ## find the best solution according to the gap method
+    function ind = gapSearch (this)
+      if (strcmp (this.SearchMethod, "globalmaxse"))
+        [gapmax, indgp] = max (this.CriterionValues);
+        for iter = 1 : length (this.InspectedK)
+          ind = iter;
+          if (this.CriterionValues(iter) > (gapmax - this.SE(indgp)))
+            return
+          endif
+        endfor
+      elseif (strcmp (this.SearchMethod, "firstmaxse"))
+        for iter = 1 : (length (this.InspectedK) - 1)
+          ind = iter;
+          if (this.CriterionValues(iter) > (this.CriterionValues(iter + 1) - ...
+                                            this.SE(iter + 1)))
+            return
+          endif
+        endfor
+      else
+        ## this should not happen
+        error (["GapEvaluation: unexpected error, please report this bug"]);
+      endif
+    endfunction
+  endmethods
+endclassdef
diff -r d354b1b49cad -r 933bc06e11e4 inst/private/SilhouetteEvaluation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/private/SilhouetteEvaluation.m	Sat Apr 03 17:51:12 2021 +0200
@@ -0,0 +1,254 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+classdef SilhouetteEvaluation < ClusterCriterion
+  ## -*- texinfo -*-
+  ## @deftypefn {Function File} {@var{eva} =} evalclusters (@var{x}, @var{clust}, @qcode{silhouette})
+  ## @deftypefnx {Function File} {@var{eva} =} evalclusters (@dots{}, @qcode{Name}, @qcode{Value})
+  ##
+  ## A silhouette object to evaluate clustering solutions.
+  ##
+  ## A @code{SilhouetteEvaluation} object is a @code{ClusterCriterion}
+  ## object used to evaluate clustering solutions using the silhouette
+  ## criterion.
+  ##
+  ## List of public properties specific to @code{SilhouetteEvaluation}:
+  ## @table @code
+  ## @item @qcode{Distance} 
+  ## a valid distance metric name, or a function handle or a numeric array as
+  ## generated by the @code{pdist} function.
+  ##
+  ## @item @qcode{ClusterPriors} 
+  ## a valid name for the evaluation of silhouette values: @code{empirical}
+  ## (default) or @code{equal}.
+  ##
+  ## @item @qcode{ClusterSilhouettes} 
+  ## a cell array with the silhouette values of each data point for each cluster
+  ## number.
+  ##
+  ## @end table
+  ## 
+  ## The best solution according to the silhouette criterion is the one
+  ## that scores the highest average silhouette value.
+  ## @end deftypefn
+  ##
+  ## @seealso{ClusterCriterion, evalclusters, silhouette}
+  
+  properties (GetAccess = public, SetAccess = private)
+    Distance = ""; # pdist parameter
+    ClusterPriors = ""; # evaluation of silhouette values: equal or empirical
+    ClusterSilhouettes = {}; # results of the silhoutte function for each K
+  endproperties
+  
+  properties (Access = protected)
+    DistanceVector = []; # vector of pdist distances
+  endproperties
+  
+  methods (Access = public)
+    ## constructor
+    function this = SilhouetteEvaluation (x, clust, KList, ...
+                    distanceMetric = "sqeuclidean", clusterPriors = "empirical")
+      this@ClusterCriterion(x, clust, KList);
+      
+      ## parsing the distance criterion
+      if (ischar (distanceMetric))
+        if (any (strcmp (lower (distanceMetric), {"sqeuclidean", ...
+                  "euclidean", "cityblock", "cosine", "correlation", ...
+                  "hamming", "jaccard"})))
+          this.Distance = lower (distanceMetric);
+          
+          ## kmeans can use only a subset
+          if (strcmp (clust, "kmeans") && any (strcmp (this.Distance, ...
+              {"euclidean", "jaccard"})))
+            error (["SilhouetteEvaluation: invalid distance criterion '%s' "...
+                    "for 'kmeans'"], distanceMetric);
+          endif
+        else
+          error ("SilhouetteEvaluation: unknown distance criterion '%s'", ...
+                 distanceMetric);
+        endif
+      elseif (isa (distanceMetric, "function_handle"))
+        this.Distance = distanceMetric;
+        
+        ## kmeans cannot use a function handle
+        if (strcmp (clust, "kmeans"))
+          error (["SilhouetteEvaluation: invalid distance criterion for "...
+                  "'kmeans'"]);
+        endif
+      elseif (isvector (distanceMetric) && isnumeric (distanceMetric))
+        this.Distance = "";
+        this.DistanceVector = distanceMetric; # the validity check is delegated
+        
+        ## kmeans cannot use a distance vector
+        if (strcmp (clust, "kmeans"))
+          error (["SilhouetteEvaluation: invalid distance criterion for "...
+                  "'kmeans'"]);
+        endif
+      else
+        error ("SilhouetteEvaluation: invalid distance metric");
+      endif
+
+      ## parsing the prior probabilities of each cluster
+      if (ischar (distanceMetric))
+        if (any (strcmp (lower (clusterPriors), {"empirical", "equal"})))
+          this.ClusterPriors = lower (clusterPriors);
+        else
+          error (["SilhouetteEvaluation: unknown prior probability criterion"...
+                 " '%s'"], clusterPriors);
+        endif
+      else
+        error ("SilhouetteEvaluation: invalid prior probabilities");
+      endif
+      
+      this.CriterionName = "silhouette";
+      this.evaluate(this.InspectedK); # evaluate the list of cluster numbers
+    endfunction
+    
+    ## set functions
+
+    ## addK
+    ## add new cluster sizes to evaluate
+    function this = addK (this, K)
+      addK@ClusterCriterion(this, K);
+      
+      ## if we have new data, we need a new evaluation
+      if (this.OptimalK == 0)
+        ClusterSilhouettes_tmp = {};
+        pS = 0; # position shift of the elements of ClusterSilhouettes 
+        for iter = 1 : length (this.InspectedK)
+          ## reorganize ClusterSilhouettes according to the new list
+          ## of cluster numbers
+          if (any (this.InspectedK(iter) == K))
+            pS += 1;
+          else
+            ClusterSilhouettes_tmp{iter} = this.ClusterSilhouettes{iter - pS};
+          endif
+        endfor
+        this.ClusterSilhouettes = ClusterSilhouettes_tmp;
+        this.evaluate(K); # evaluate just the new cluster numbers
+      endif
+    endfunction
+    
+    ## compact
+    ## ...
+    function this = compact (this)
+      # FIXME: stub!
+      warning ("SilhouetteEvaluation: compact is unavailable");
+    endfunction
+  endmethods
+  
+  methods (Access = protected)
+    ## evaluate
+    ## do the evaluation
+    function this = evaluate (this, K)
+      ## use complete observations only
+      UsableX = this.X(find (this.Missing == false), :);
+      if (! isempty (this.ClusteringFunction))
+        ## build the clusters
+        for iter = 1 : length (this.InspectedK)
+          ## do it only for the specified K values
+          if (any (this.InspectedK(iter) == K))
+            if (isa (this.ClusteringFunction, "function_handle"))
+              ## custom function
+              ClusteringSolution = ...
+                this.ClusteringFunction(UsableX, this.InspectedK(iter));
+              if (ismatrix (ClusteringSolution) && ...
+                  rows (ClusteringSolution) == this.NumObservations && ...
+                  columns (ClusteringSolution) == this.P)
+                ## the custom function returned a matrix:
+                ## we take the index of the maximum value for every row
+                [~, this.ClusteringSolutions(:, iter)] = ...
+                  max (ClusteringSolution, [], 2);
+              elseif (iscolumn (ClusteringSolution) &&
+                      length (ClusteringSolution) == this.NumObservations)
+                this.ClusteringSolutions(:, iter) = ClusteringSolution;
+              elseif (isrow (ClusteringSolution) &&
+                      length (ClusteringSolution) == this.NumObservations)
+                this.ClusteringSolutions(:, iter) = ClusteringSolution';
+              else
+                error (["SilhouetteEvaluation: invalid return value from " ...
+                        "custom clustering function"]);
+              endif
+              this.ClusteringSolutions(:, iter) = ...
+                this.ClusteringFunction(UsableX, this.InspectedK(iter));
+            else
+              switch (this.ClusteringFunction)
+                case "kmeans"
+                  this.ClusteringSolutions(:, iter) = kmeans (UsableX, ...
+                    this.InspectedK(iter), "Distance", this.Distance, ...
+                    "EmptyAction", "singleton", "Replicates", 5);
+                case "linkage"
+                  if (! isempty (this.Distance))
+                    ## use clusterdata
+                    Distance_tmp = this.Distance;
+                    LinkageMethod = "average"; # for non euclidean methods
+                    if (strcmp (this.Distance, "sqeuclidean"))
+                      ## pdist uses different names for its algorithms
+                      Distance_tmp = "squaredeuclidean";
+                      LinkageMethod = "ward";
+                    elseif (strcmp (this.Distance, "euclidean"))
+                      LinkageMethod = "ward";
+                    endif
+                    this.ClusteringSolutions(:, iter) = clusterdata (UsableX,...
+                      "MaxClust", this.InspectedK(iter), ...
+                      "Distance", Distance_tmp, "Linkage", LinkageMethod);
+                  else
+                    ## use linkage
+                    Z = linkage (this.DistanceVector, "average");
+                    this.ClusteringSolutions(:, iter) = ...
+                         cluster (Z, "MaxClust", this.InspectedK(iter));
+                  endif
+                case "gmdistribution"
+                  gmm = fitgmdist (UsableX, this.InspectedK(iter), ...
+                        "SharedCov", true, "Replicates", 5);
+                  this.ClusteringSolutions(:, iter) = cluster (gmm, UsableX);
+                otherwise
+                  error (["SilhouetteEvaluation: unexpected error, " ...
+                         "report this bug"]);
+              endswitch
+            endif
+          endif
+        endfor
+      endif
+      
+      ## get the silhouette values for every clustering
+      set (0, 'DefaultFigureVisible', 'off'); # temporarily disable figures
+      for iter = 1 : length (this.InspectedK)
+        ## do it only for the specified K values
+        if (any (this.InspectedK(iter) == K))
+          this.ClusterSilhouettes{iter} = silhouette (UsableX, ...
+                                            this.ClusteringSolutions(:, iter));
+          if (strcmp (this.ClusterPriors, "empirical"))
+            this.CriterionValues(iter) = mean (this.ClusterSilhouettes{iter});
+          else
+            ## equal
+            this.CriterionValues(iter) = 0;
+            si = this.ClusterSilhouettes{iter};
+            for k = 1 : this.InspectedK(iter)
+              this.CriterionValues(iter) += mean (si(find ...
+                                     (this.ClusteringSolutions(:, iter) == k))); 
+            endfor
+            this.CriterionValues(iter) /= this.InspectedK(iter);
+          endif
+        endif
+      endfor
+      set (0, 'DefaultFigureVisible', 'on'); # enable figures again
+      
+      [~, this.OptimalIndex] = max (this.CriterionValues);
+      this.OptimalK = this.InspectedK(this.OptimalIndex(1));
+      this.OptimalY = this.ClusteringSolutions(:, this.OptimalIndex(1));
+    endfunction
+  endmethods
+endclassdef
