# HG changeset patch
# User Stefano Guidoni <ilguido@users.sf.net>
# Date 1614623252 -3600
#      Mon Mar 01 19:27:32 2021 +0100
# Node ID 041359224c4e4228d71e41dab4d4118bc2af493d
# Parent  c7da8294729ff5f53b975a60cf37900bf834d13f
update function dendogram

diff -r c7da8294729f -r 041359224c4e inst/dendrogram.m
--- a/inst/dendrogram.m	Sat Feb 27 18:21:40 2021 +0100
+++ b/inst/dendrogram.m	Mon Mar 01 19:27:32 2021 +0100
@@ -1,109 +1,334 @@
-%% Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%%    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
-%%    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/>.
+## Copyright (C) 2021 Stefano Guidoni <ilguido@users.sf.net>
+## Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+##
+## 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{p} = } dendrogram (@var{tree})
-%% @deftypefnx {Function File} {@var{p, t} = } dendrogram (@var{tree})
-%% @deftypefnx {Function File} {@var{p, t, perm} = } dendrogram (@var{tree})
-%% Plots a dendrogram using the output of function @command{linkage}.
-%%
-%% t is a vector containing the leaf node number for each object in the
-%% original dataset.  For now, all objects are leaf nodes.
-%%
-%% perm is the permutation of the input objects used to display the
-%% dendrogram, in left-to-right order.
-%%
-%% TODO: Return handle to lines to set properties
-%% TODO: Rescale the plot automatically based on data.
-%%
-%% @seealso{linkage}
-%% @end deftypefn
+## -*- texinfo -*-
+## @deftypefn {Function File} {} dendrogram (@var{tree})
+## @deftypefnx {Function File} {} dendrogram (@var{tree}, @var{p})
+## @deftypefnx {Function File} {} dendrogram (@var{tree}, @var{prop}, @var{val})
+## @deftypefnx {Function File} @
+##   {} dendrogram (@var{tree}, @var{p}, @var{prop}, @var{val} )
+## @deftypefnx {Function File} {@var{h} = } dendrogram (@dots{})
+## @deftypefnx {Function File} @
+##   {[@var{h}, @var{t}, @var{perm}] = } dendrogram (@dots{})
+##
+## Plot a dendrogram of a hierarchical binary cluster tree.
+##
+## Given @var{tree}, a hierarchical binary cluster tree as the output of
+## @code{linkage}, plot a dendrogram of the tree. The number of leaves shown by
+## the dendrogram plot is limited to @var{p}. The default value for @var{p} is 
+## 30. Set @var{p} to 0 to plot all leaves.
+##
+## The optional outputs are @var{h}, @var{t} and @var{perm}: 
+## @itemize @bullet
+## @item @var{h} is a handle to the lines of the plot.
+##
+## @item @var{t} is the vector with the numbers assigned to each leaf. 
+## Each element of @var{t} is a leaf of @var{tree} and its value is the number 
+## shown in the plot.
+## When the dendrogram plot is collapsed, that is when the number of shown
+## leaves @var{p} is inferior to the total number of leaves, a single leaf of
+## the plot can represent more than one leaf of @var{tree}: in that case
+## multiple elements of @var{t} share the same value, that is the same leaf of
+## the plot.
+## When the dendrogram plot is not collapsed, each leaf of the plot is the leaf
+## of @var{tree} with the same number.
+##
+## @item @var{perm} is the vector list of the leaves as ordered as in the plot.
+## @end itemize
+##
+## Additional input properties can be specified by pairs of properties and 
+## values. Known properties are:
+## @itemize @bullet
+## @item @qcode{"Reorder"}
+## Reorder the leaves of the dendrogram plot using a numerical vector of size n,
+## the number of leaves. When @var{p} is smaller than @var{n}, the reordering
+## cannot break the @var{p} groups of leaves.
+##
+## @item @qcode{"Orientation"}
+## Change the orientation of the plot. Available values: @qcode{top} (default), 
+## @qcode{bottom}, @qcode{left}, @qcode{right}.
+##
+## @item @qcode{"CheckCrossing"}
+## Not implemented.
+##
+## @item @qcode{"ColorThreshold"}
+## Not implemented.
+##
+## @item @qcode{"Labels"}
+## Not implemented.
+## @end itemize
+##
+## @end deftypefn
+##
+## @seealso{cluster, clusterdata, cophenet, inconsistent,linkage, pdist}
 
-function [p, t, perm] = dendrogram (tree)
+function [H, T, perm] = dendrogram (tree, varargin)
 
   [m d] = size (tree);
-  if d != 3
-    error ("Input data must be a tree as returned by function linkage.")
+  if ((d != 3) || (! isnumeric (tree)) || ...
+      (! (max (tree(end,1:2)) == m * 2)))
+    error (["dendrogram: tree must be a matrix as generated by the " ...
+      "linkage function"]);
   end
-  n = m + 1;
-
-  % t is the leaf node number for all objects in the original dataset.
-  % TODO: Add support for collapsing the tree.
-  % For now, we always display all objects, so this is the identity map.
-  t = (1:m)';
-
-  nc = max(tree(:,1:2)(:));
+  
+  pair_index = 1;
 
-  % Vector with the horizontal and vertical position of each cluster
-  p = zeros (nc,2);
-
-  perm = zeros (n,1);
+  ## node count
+  n = m + 1;
+    
+  P = 30; # default value
+  vReorder = []; 
+  orientation = "top";
+  if (nargin > 1)
+    if (isnumeric (varargin{1}) && isscalar (varargin{1})) 
+      ## dendrogram (tree, P)
+      P = varargin{1};
+      pair_index++;
+    endif
+    
+    ## dendrogram (..., Name, Value)
+    while (pair_index < (nargin - 1))
+      switch (lower (varargin{pair_index}))
+        case "reorder"
+          if (isvector (varargin{pair_index + 1}) && ...
+            isnumeric (varargin{pair_index + 1}) && ...
+            length (varargin{pair_index + 1}) == n )
+            vReorder = varargin{pair_index + 1};
+          else
+            error (["dendrogram: 'reorder' must be a numeric vector of size" ...
+              "n, the number of leaves"]);
+          endif
+        case "checkcrossing"
+          warning ("dendrogram: property '%s' not implemented",...
+            varargin{pair_index});
+        case "colorthreshold"
+          warning ("dendrogram: property '%s' not implemented",...
+            varargin{pair_index});
+        case "orientation"
+          orientation = varargin{pair_index + 1}; # validity check below
+        case "labels"
+          warning ("dendrogram: property '%s' not implemented",...
+           varargin{pair_index});
+        otherwise
+          error ("dendrogram: unknown property '%s'", varargin{pair_index});
+      endswitch
+      
+      pair_index += 2;
+    endwhile 
+  endif
 
-  %% Ordering by depth-first search
-  nodecount = 0;
-  nodes_to_visit = nc+1;
-  while !isempty(nodes_to_visit)
-    currentnode = nodes_to_visit(1);
-    nodes_to_visit(1) = [];
-    if currentnode > n
-      node = currentnode - n;
-      nodes_to_visit = [tree(node,[2 1]) nodes_to_visit];
-    end
+  ## MATLAB compatibility:
+  ## P <= 0 to plot all leaves
+  if P < 1
+    P = n;
+  endif
+  
+  if n > P
+    level_0 = tree((n - P),3);
+  else
+    P = n;
+    level_0 = 0;
+  endif
 
-    if currentnode <= n && p(currentnode,1) == 0
-      nodecount +=1;
-      p(currentnode,1) = nodecount;
-      perm(nodecount) = currentnode;
-    end
-
-  end
-
-  % Compute the horizontal position, begin-end
-  % and vertical position of all clusters.
-  for i = 1:m
-    p(n+i,1)   = mean (p(tree(i,1:2),1));
-    p(n+i,2)   = tree(i,3);
-    x(i,1:2) = p(tree(i,1:2),1);
-  end
-
+  vLeafPosition = zeros((n + m), 1);
+  T = (1:n)';
+  nodecnt = 1;
+  
+  ## main
+  dendrogram_recursive (m, 0);
+  
+  ## T reordering
+  ## MATLAB compatibility: when n > P, each node group is renamed with a number
+  ## between 1 and P, according to the smallest node index of each group;
+  ## the group with the node 1 is always group 1, while group 2 is the group
+  ## with the smallest node index outside of group 1, and group 3 is the group
+  ## with the smallest node index outside of groups 1 and 2...
+  newT = 1 : (length (T));
+  if n > P
+    uniqueT = unique (T);
+    minT = zeros (uniqueT, 1);
+    counter = 1;
+    
+    for i = 1 : length (uniqueT) # it should be exactly equal to P
+      idcs = find (T == uniqueT(i));
+      minT(i) = min (idcs);
+    endfor
+    
+    minT = minT(find (minT > 0)); # to prevent a strange bug
+    [minT, minTidcs] = sort (minT);
+    uniqueT = uniqueT(minTidcs);
+    for i = 1 : length (uniqueT)
+      idcs = find (T == uniqueT(i));
+      newT(idcs) = counter++;
+    endfor
+  endif
+  
+  ## leaf reordering
+  if (! isempty(vReorder))
+    if (P < n)
+      checkT = newT(vReorder(:));
+      for i = 1 : P
+        idcs = find (checkT == i);
+        if (length (idcs) > 1)
+          if (max (idcs) - min (idcs) >= length (idcs))
+            error (["dendrogram: invalid reordering that redefines the 'P'"...
+              "groups of leaves"]);
+          endif
+        endif
+      endfor
+      checkT = unique (checkT, "stable");
+      vNewLeafPosition = zeros (n, 1);
+      uT = unique (T, "stable");
+      for i = 1 : P
+        vNewLeafPosition(uT(checkT(i))) = i;
+      endfor
+      vLeafPosition = vNewLeafPosition;
+    else
+      for i = 1 : length (vReorder)
+        vLeafPosition(vReorder(i)) = i;
+      endfor
+    endif
+  endif
+    
+  ## figure
+  x = [];
+  
   figure(gcf)
-  % plot horizontal lines
-  tmp = line (x', tree(:,[3 3])');
-
-  % plot vertical lines
-  [~,tf]  = ismember (1:nc, tree(:,1:2));
-  [ind,~] = ind2sub (size (tree(:,1:2)), tf);
-  y       = [p(1:nc,2) tree(ind,3)];
-  tmp     = line ([p(1:nc,1) p(1:nc,1)]',y');
-
-  xticks = 1:n;
-  xl_txt = arrayfun (@num2str, perm,"uniformoutput",false);
-  set (gca,"xticklabel",xl_txt,"xtick",xticks);
-  axis ([0.5 n+0.5 0 max(tree(:,3))+0.1*min(tree(:,3))]);
-
+  
+  ## ticks and tricks
+  xticks = 1:P;
+  perm = zeros (P, 1);
+  for i = 1 : length (vLeafPosition)
+    if (vLeafPosition(i) != 0)
+        idcs = find (T == i);
+        perm(vLeafPosition(i)) = newT(idcs(1));
+    endif
+  endfor
+  T = newT; # this should be unnecessary for n <= P
+  
+  ## lines
+  for i = (n - P + 1) : m
+    vLeafPosition(n + i) = mean (vLeafPosition(tree(i,1:2),1));
+    x(end + 1,1:4) = [vLeafPosition(tree(i,1:2))' tree(i,[3 3])];
+    for j = 1 : 2
+      x0 = 0;
+      if tree(i,j) > (2 * n - P)
+        x0 = tree(tree(i,j) - n,3);
+      endif
+      x(end + 1,1:4) = [vLeafPosition(tree(i,[j j]))' x0 tree(i,3)];  
+    endfor    
+  endfor
+  
+  ## plot stuff
+  if (strcmp (orientation, "top"))
+    H = line (x(:,1:2)', x(:,3:4)', "color", "blue");
+    set (gca,"xticklabel",perm,"xtick",xticks);
+  elseif (strcmp (orientation, "bottom"))
+    H = line (x(:,1:2)', x(:,3:4)', "color", "blue");
+    set (gca,"xticklabel",perm,"xtick",xticks,"xaxislocation", "top");
+    axis ("ij");
+  elseif (strcmp (orientation, "left"))
+    H = line (x(:,3:4)', x(:,1:2)', "color", "blue");
+    set (gca,"yticklabel",perm,"ytick",xticks,"xdir","reverse",...
+         "yaxislocation", "right");
+  elseif (strcmp (orientation, "right"))
+    H = line (x(:,3:4)', x(:,1:2)', "color", "blue");
+    set (gca,"yticklabel",perm,"ytick",xticks);
+  else
+    error ("dendrogram: invalid orientation '%s'", orientation);
+  endif
+  
+  ## dendrogram_recursive
+  function dendrogram_recursive (k, cn)
+    if (tree(k,3) > level_0)
+      for j = 1:2
+        if (tree(k,j) > n)
+          dendrogram_recursive (tree(k,j) - n, 0)
+        else
+          vLeafPosition(tree(k,j)) = nodecnt++;
+          T(tree(k,j)) = tree(k,j);
+        endif
+      endfor
+    else
+      for j = 1:2
+        if (cn == 0)
+          cn = n + k;
+          vLeafPosition(cn) = nodecnt++;
+        endif
+        
+        if (tree(k,j) > n)
+          dendrogram_recursive (tree(k,j) - n, cn)
+        else
+          T(tree(k,j)) = cn;
+        endif
+      endfor
+    endif
+  endfunction
 endfunction
 
-%!demo
+## Test input validation
+%!error dendrogram ()
+%!error <tree must be .*> dendrogram (ones (2,2), 1)
+%!error <unknown property .*> dendrogram ([1 2 1], 1, "xxx", "xxx")
+%!error <reorder.*> dendrogram ([1 2 1], "Reorder", "xxx")
+%!error <reorder.*> dendrogram ([1 2 1], "Reorder", [1 2 3 4])
+%!error <invalid orientation .*> dendrogram ([1 2 1], "Orientation", "north")
+## Demonstrations
+## 1.
+%!demo 1
 %! y      = [4 5; 2 6; 3 7; 8 9; 1 10];
 %! y(:,3) = 1:5;
-%! figure(gcf); clf;
-%! dendrogram(y);
-
-%!demo
-%! v = 2*rand(30,1)-1;
-%! d = abs(bsxfun(@minus, v(:,1), v(:,1)'));
+%! figure (gcf); clf;
+%! dendrogram (y);
+## 2.
+%!demo 2
+%! v = 2 * rand (30,1) - 1;
+%! d = abs (bsxfun (@minus, v(:,1), v(:,1)'));
 %! y = linkage (squareform(d,"tovector"));
-%! figure(gcf); clf;
-%! dendrogram(y);
+%! figure (gcf); clf;
+%! dendrogram (y);
+## 3. collapsed tree
+%!demo "collapsed tree, find all the leaves of node 5"
+%! X = randn (60,2);
+%! D = pdist (X);
+%! y = linkage (D, "average");
+%! figure (gcf); clf;
+%! subplot (2, 1, 1);
+%! title ("original tree");
+%! dendrogram (y, 0);
+%! subplot (2, 1, 2);
+%! title ("collapsed tree");
+%! [~,t] = dendrogram (y, 20);
+%! find(t == 5)
+## 4. optimal leaf order
+%!demo "optimal leaf order"
+%! X = randn (30,2);
+%! D = pdist (X);
+%! y = linkage (D, "average");
+%! order = optimalleaforder (y, D);
+%! figure (gcf); clf;
+%! subplot (2, 1, 1);
+%! title ("original leaf order");
+%! dendrogram (y);
+%! subplot (2, 1, 2);
+%! title ("optimal leaf order");
+%! dendrogram (y, "Reorder", order);
+## 5. orientation
+%!demo "horizontal orientation"
+%! X = randn (10,2);
+%! D = pdist (X);
+%! y = linkage (D, "average");
+%! dendrogram (y, "Orientation", "left");
