# HG changeset patch
# User Stefano Guidoni <ilguido@users.sf.net>
# Date 1615462590 -3600
#      Thu Mar 11 12:36:30 2021 +0100
# Node ID 76bdc63f455eb774f8aa0d00f18fba06d3011b3a
# Parent  b5a88871519cbad1a2bd2dac4dfff75f351cd5ee
silhouette: new function

diff -r b5a88871519c -r 76bdc63f455e INDEX.in
--- a/INDEX.in	Sun Mar 07 21:29:34 2021 +0100
+++ b/INDEX.in	Thu Mar 11 12:36:30 2021 +0100
@@ -79,6 +79,7 @@
  hist3
  normplot
  repanova
+ silhouette
  violin
 Models
  hmmestimate 
diff -r b5a88871519c -r 76bdc63f455e NEWS
--- a/NEWS	Sun Mar 07 21:29:34 2021 +0100
+++ b/NEWS	Thu Mar 11 12:36:30 2021 +0100
@@ -20,6 +20,8 @@
  
  ** optimalleaforder.m (patch #10034)
  
+ ** silhouette.m       (patch #9743)
+ 
  ** slicesample        (patch #10019)
 
  Improvements:
diff -r b5a88871519c -r 76bdc63f455e inst/silhouette.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/silhouette.m	Thu Mar 11 12:36:30 2021 +0100
@@ -0,0 +1,200 @@
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+## Copyright (C) 2016 Nan Zhou <zhnanx@gmail.com>
+##
+## 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{si}, @var{h}] =} silhouette (@var{X}, @var{clust})
+## @deftypefnx {Function File} @
+##  {[@var{si}, @var{h}] =} silhouette (@dots{}, @var{Metric}, @var{MetricArg})
+##
+## Compute the silhouette values of clustered data and show them on a plot.
+##
+## @var{X} is a n-by-p matrix of n data points in a p-dimensional space.  Each
+## datapoint is assigned to a cluster using @var{clust}, a vector of n elements,
+## one cluster assignment for each data point.
+## 
+## Each silhouette value of @var{si}, a vector of size n, is a measure of the
+## likelihood that a data point is accurately classified to the right cluster.  
+## Defining "a" as the mean distance between a point and the other points from 
+## its cluster, and "b" as the mean distance between that point and the points 
+## from other clusters, the silhouette value of the i-th point is:
+##
+## @tex
+## \def\frac#1#2{{\begingroup#1\endgroup\over#2}}
+## $$ S_i = \frac{b_i - a_i}{max(a_1,b_i)} $$
+## @end tex
+## @ifnottex
+## @verbatim
+##          bi - ai  
+## Si =  ------------
+##        max(ai,bi)    
+## @end verbatim
+## @end ifnottex
+##
+## Each element of @var{si} ranges from -1, minimum likelihood of a correct
+## classification, to 1, maximum likelihood.
+##
+## Optional input value @var{Metric} is the metric used to compute the distances
+## between data points. Since @code{silhouette} uses @code{pdist} to compute
+## these distances, @var{Metric} is quite similar to the option @var{Metric} of
+## pdist and it can be:
+## @itemize @bullet
+## @item A known distance metric defined as a string: @qcode{Euclidean},
+## @qcode{sqEuclidean} (default), @qcode{cityblock}, @qcode{cosine},
+## @qcode{correlation}, @qcode{Hamming}, @qcode{Jaccard}.
+##
+## @item A vector as those created by @code{pdist}. In this case @var{X} does
+## nothing.
+##
+## @item A function handle that is passed to @code{pdist} with @var{MetricArg}
+## as optional inputs.
+## @end itemize
+##
+## Optional return value @var{h} is a handle to the silhouette plot.
+##
+## @strong{Reference}
+## Peter J. Rousseeuw, Silhouettes: a Graphical Aid to the Interpretation and 
+## Validation of Cluster Analysis. 1987. doi:10.1016/0377-0427(87)90125-7
+## @end deftypefn
+##
+## @seealso{dendrogram, evalcluster, kmeans, linkage, pdist}
+
+function [si, h] = silhouette (X, clust, metric = "sqeuclidean", varargin)
+  ## check the input parameters
+  if (nargin < 2)
+    print_usage ();
+  endif
+  
+  n = size (clust, 1);
+  
+  ## check size
+  if (! isempty (X))
+    if (size (X, 1) != n)
+      error ("First dimension of X <%d> doesn't match that of clust <%d>",...
+        size (X, 1), n);
+    endif
+  endif
+  
+  ## check metric
+  if (ischar (metric))
+    metric = lower (metric);
+    switch (metric)
+      case "sqeuclidean"
+        metric = "seuclidean";
+      case { "euclidean", "cityblock", "cosine", ...
+           "correlation", "hamming", "jaccard" }
+        ;
+      otherwise
+        error ("silhouette: invalid metric '%s'", metric);
+    endswitch
+  elseif (isnumeric (metric) && isvector (metric))
+    ## X can be omitted when using this
+    distMatrix = squareform (metric);
+    if (size (distMatrix, 1) != n)
+      error ("First dimension of X <%d> doesn't match that of clust <%d>",...
+        size (distMatrix, 1), n);
+    endif
+  endif
+  
+  ## main
+  si = zeros(n, 1);
+  clusterIDs = unique (clust); # eg [1; 2; 3; 4]
+  m = length (clusterIDs);
+
+  ## distance matrix showing the distance for any two rows of X
+  if (! exist ('distMatrix', 'var'))
+    distMatrix = squareform (pdist (X, metric, varargin));
+  endif
+  
+  ## calculate values of si one by one
+  for iii = 1 : length (si)
+    
+    ## distance of iii to all others
+    iii2all = distMatrix(iii, :);
+    
+    ## allocate values to clusters
+    groupedValues = {};
+    for jjj = 1 : m
+      groupedValues{clusterIDs(jjj)} = [iii2all(clust == clusterIDs(jjj))];
+    endfor
+    ## end allocation
+    
+    ## calculate a(i)
+    ## dist of object iii to all other objects in the same cluster
+    a_iii = groupedValues{clust(iii)};
+    ## average distance of iii to all other objects in the same cluster
+    if (length (a_iii) == 1)
+      si(iii) = 0;
+      continue;
+    else
+      a_i = (sum (a_iii)) / (length (a_iii) - 1);
+    endif
+    ## end a(i)
+    
+    
+    ## calculate b(i)
+    clusterIDs_new = clusterIDs;
+    ## remove the cluster iii in
+    clusterIDs_new(find (clusterIDs_new == clust(iii))) = [];
+    ## average distance of iii to all objects of another cluster
+    a_iii_2others = zeros (length (clusterIDs_new), 1);
+    for jjj = 1 : length (clusterIDs_new)
+      values_another = groupedValues{clusterIDs_new(jjj)};
+      a_iii_2others(jjj) = mean (values_another);
+    endfor
+    b_i = min (a_iii_2others);
+    ## end b(i)
+    
+    
+    ## calculate s(i)
+    si(iii) = (b_i - a_i) / (max ([a_i; b_i]));
+    ## end s(i)
+  endfor
+
+  ## plot
+  ## a poor man silhouette graph
+  vBarsc = zeros (m, 1);
+  vPadding = [0; 0; 0; 0];
+  Bars = vPadding;
+  
+  for i = 1 : m
+    vBar = si(find (clust == clusterIDs(i)));
+    vBarsc(i) = length (Bars) + (length (vBar) / 2);
+    Bars = [Bars; (sort (vBar, "descend")); vPadding];
+  endfor
+  
+  h = barh (Bars, "hist", "facecolor", [0 0.4471 0.7412]);
+  
+  xlabel ("Silhouette Value");
+  ylabel ("Cluster");
+  set (gca, "ytick", vBarsc, "yticklabel", clusterIDs);
+  ylim ([0 (length (Bars))]);
+  axis ("ij");
+endfunction
+ 
+%!error silhouette ();
+%!error silhouette ([1 2; 1 1]);
+%!error <X .* doesn't match .* clust> silhouette ([1 2; 1 1], [1 2 3]');
+%!error <invalid metric> silhouette ([1 2; 1 1], [1 2]', "xxx");
+
+%!demo
+%! load fisheriris;
+%! X = meas(:,3:4);
+%! cidcs = kmeans (X, 3, "Replicates", 5);
+%! silhouette (X, cidcs);
+%! y_labels(cidcs([1 51 101])) = unique (species);
+%! set (gca, "yticklabel", y_labels);
+%! title ("Fisher's iris data");
