## Copyright (C) 2021 Stefano Guidoni ## ## 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; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {} {@var{leafOrder} =} optimalleaforder (@var{tree}, @var{D}) ## @deftypefnx {Function File} @ ## {@var{leafOrder} =} optimalleaforder (..., @var{Name}, @var{Value}) ## ## Compute the optimal leaf ordering of a hierarchical binary cluster tree. ## ## The optimal leaf ordering of a tree is the ordering which minimizes the sum ## of the distances between each leaf and its adjacent leaves, without altering ## the structure of the tree, that is without redefining the clusters of the ## tree. ## ## Required inputs: ## @itemize ## @item ## @var{tree}: a hierarchical cluster tree @var{tree} generated by the ## @code{linkage} function. ## ## @item ## @var{D}: a matrix of distances as computed by @code{pdist}. ## @end itemize ## ## Optional inputs can be the following property/value pairs: ## @itemize ## @item ## property 'Criteria' at the moment can only have the value 'adjacent', ## for minimizing the distances between leaves. ## ## @item ## property 'Transformation' can have one of the values 'linear', 'inverse' ## or a handle to a custom function which computes @var{S} the similarity ## matrix. ## @end itemize ## ## optimalleaforder's output @var{leafOrder} is the optimal leaf ordering. ## ## @strong{Reference} ## Bar-Joseph, Z., Gifford, D.K., and Jaakkola, T.S. Fast optimal leaf ordering ## for hierarchical clustering. Bioinformatics vol. 17 suppl. 1, 2001. ## @end deftypefn ## ## @seealso{dendrogram,linkage,pdist} function leafOrder = optimalleaforder ( varargin ) ## check the input if ( nargin < 2 ) print_usage (); endif tree = varargin{1}; D = varargin{2}; criterion = "adjacent"; # default and only value at the moment transformation = "inverse"; if ((columns (tree) != 3) || (! isnumeric (tree)) || ... (! (max (tree(end, 1:2)) == rows (tree) * 2))) error (["optimalleaforder: tree must be a matrix as generated by the " ... "linkage function"]); endif ## read the paired arguments if (! all (cellfun ("ischar", varargin(3:end)))) error ("optimalleaforder: character inputs expected for arguments 3 and up"); else varargin(3:end) = lower (varargin(3:end)); endif pair_index = 3; while (pair_index <= (nargin - 1)) switch (varargin{pair_index}) case "criteria" criterion = varargin{pair_index + 1}; if (strcmp (criterion, "group")) ## MATLAB compatibility: ## the 'group' criterion is not implemented error ("optimalleaforder: unavailable criterion 'group'"); elseif ! strcmp (criterion, "adjacent") error ("optimalleaforder: invalid criterion %s", criterion); endif case "transformation" transformation = varargin{pair_index + 1}; otherwise error ("optimalleaforder: unknown property %s", varargin{pair_index}); endswitch pair_index += 2; endwhile ## D can be either a vector or a matrix, ## but it is easier to work with a matrix if (isvector (D)) D = squareform (D); endif n = rows (D); m = rows (tree); if (n != (m + 1)) error (["optimalleaforder: D must be a matrix or vector generated by " ... "the pdist function"]); endif ## the similarity matrix, basically an inverted distance matrix S = zeros (n); if (strcmpi (transformation, "linear")) ## linear similarity maxD = max (max (D)); S = maxD .- D; elseif (strcmpi (transformation, "inverse")) ## similarity as inverted distance S = 1 ./ D; elseif (is_function_handle (transformation)) ## custom similarity S = feval (transformation, D); else error ("optimalleaforder: invalid transformation %s", transformation); endif ## main body ## for each node v we compute the maximum similarity of the subtree M(w,u,v), ## where the leftmost leaf is w and the rightmost is u; remember that ## M(w,u,v) = M(u,w,v) M = zeros (n, n, n + m); ## O is a utility matrix: for each node of the tree we store the left and ## right leaves of the optimal subtree O = [1:( n + m ); 1:( n + m ); (zeros (1, (n + m)))]'; ## compute M for every node v for iter = 1 : m v = iter + n; # current node l = optimalleaforder_getLeafList (tree(iter, 1)); # the left subtree r = optimalleaforder_getLeafList (tree(iter, 2)); # the right subtree if (tree(iter,1) > n) l_l = optimalleaforder_getLeafList (tree(tree(iter, 1) - n, 1)); l_r = optimalleaforder_getLeafList (tree(tree(iter, 1) - n, 2)); else l_l = l_r = l; endif if (tree(iter,2) > n) r_l = optimalleaforder_getLeafList (tree(tree(iter, 2) - n, 1)); r_r = optimalleaforder_getLeafList (tree(tree(iter, 2) - n, 2)); else r_l = r_r = r; endif ## let's find the maximum value of M(w,u,v) when: w is a leaf of the left ## subtree of v and u is a leaf of the right subtree of v for i = 1 : length (l) if (isempty (find (l(i) == l_l))) x = l_l; else x = l_r; endif for j = 1 : length (r) if (isempty (find (r(j) == r_l))) y = r_l; else y = r_r; endif ## max(M(w,u,v)) = max(M(w,k,v_l)) + max(M(h,u,v_r)) + S(k,h) ## where: v_l is the left child of v and v_r the right child of v M_tmp = repmat (M(l(i), x(:), tree(iter, 1)), length (y), 1) + ... repmat (M(y(:), r(j), tree(iter, 2)), 1, length (x)) + ... S(y(:), x(:)); M_max = max (max (M_tmp)); # this is M(l(i), r(j), v) [h, k] = find (M_tmp == M_max); M(l(i), r(j), v) = M_max; M(r(j), l(i), v) = M(l(i), r(j), v); if (M_max > O(v,3)) O(v, 1) = l(i); # this is w O(v, 2) = r(j); # this is u O(v, 3) = M_max; # this is M(w, u, v) endif endfor endfor endfor ## reordering: ## we found the M(w,u,v) corresponding to the optimal leaf order, now we can ## compute the optimal leaf order given our M(w,u,v) ## the return value leafOrder = zeros ( 1, n ); leafOrder(1) = O(end, 1); leafOrder(n) = O(end, 2); ## the inverse operation, only easier, to get the leaf order: now we know the ## leftmost and rightmost leaves of the best subtree, we may have to flip it ## though for iter = m : -1 : 1 v = iter + n; extremes = O(v, [1, 2]); l_node = tree(iter, 1); r_node = tree(iter, 2); l = optimalleaforder_getLeafList (l_node); r = optimalleaforder_getLeafList (r_node); if (l_node > n) l_l = optimalleaforder_getLeafList (tree(l_node - n, 1)); l_r = optimalleaforder_getLeafList (tree(l_node - n, 2)); else l_l = l_r = l; endif if (r_node > n) r_l = optimalleaforder_getLeafList (tree(r_node - n, 1)); r_r = optimalleaforder_getLeafList (tree(r_node - n, 2)); else r_l = r_r = r; endif ## this means that we need to flip the subtree if (isempty (find (extremes(1) == l))) l_tmp = l; l_l_tmp = l_l; l_r_tmp = l_r; l = r; l_l = r_l; l_r = r_r; r = l_tmp; r_l = l_l_tmp; r_r = l_r_tmp; node_tmp = l_node; l_node = r_node; r_node = node_tmp; endif if (isempty (find (extremes(1) == l_l))) x = l_l; else x = l_r; endif if (isempty (find (extremes(2) == r_l))) y = r_l; else y = r_r; endif M_tmp = repmat (M(extremes(1), x(:), l_node), length (y), 1) + ... repmat (M(y(:), extremes(2), r_node), 1, length (x)) + ... S(y(:), x(:)); M_max = max (max (M_tmp)); [h, k] = find (M_tmp == M_max); O(l_node, 1) = extremes(1); O(l_node, 2) = x(k); O(r_node, 1) = y(h); O(r_node, 2) = extremes(2); p_1 = find (leafOrder == extremes(1)); p_2 = find (leafOrder == extremes(2)); leafOrder (p_1 + (length (l)) - 1) = x(k); leafOrder (p_1 + (length (l))) = y(h); endfor ## function: optimalleaforder_getLeafList ## get the list of leaves under a given node function vector = optimalleaforder_getLeafList (nodes_to_visit) vector = []; 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]; endif if (currentnode <= n) vector = [vector currentnode]; endif endwhile endfunction endfunction ## Test input validation %!error optimalleaforder () %!error optimalleaforder (1) %!error optimalleaforder (ones (2, 2), 1) %!error optimalleaforder ([1 2 3], [1 2; 3 4], "criteria", 5) %!error optimalleaforder ([1 2 1], [1 2 3]) %!error optimalleaforder ([1 2 1], 1, "xxx", "xxx") %!error optimalleaforder ([1 2 1], 1, "Transformation", "xxx") ## Demonstration %!demo %! X = randn (10, 2); %! D = pdist (X); %! tree = linkage(D, 'average'); %! optimalleaforder (tree, D, 'Transformation', 'linear')