function [h,T] = dendrogram(Z,D,p)
%DENDROGRAM Generate dendragram plot for temporal sequence data.
% The Matlab 5.3 'dendrogram' routine was modified by Uri Alon,
% Oct 24, 2000 to allow temporal ordering of the dendrogram.
% The input Z is the output of the matlab routine LINKAGE.
% D is a matrix containing the log of the data. Before
% taking the log, every time dependent data sequence should be
% normalized by its maximal value. Time points are assumed to
% be equally spaced for all measurements.
% DENDROGRAM generates a dendrogram from the output matrix of
% LINKAGE, Z. Z is a (M-1) by 3 matrix. M is the number of
% observations in the original data.
%
% A dendrogram consists of many upsidedown U shape lines connecting
% nodes in a hierachichale tree. Except for the WARD linkage (see
% LINKAGE), the height of each U is the distance between the two
% clusters to be connected at that time.
%
% DENDROGRAM(Z,P) generate a dendrogram with only the top P nodes.
% When there are more than 30 nodes in the original data, the
% dendrogram may look crowded. The default value of P is 30. If P =
% 0, then, every node will be displayed.
%
% H = DENDROGRAM(...) returns a vector of line handles.
%
% [H, T] = DENDROGRAM(...) also returns T, a vector of size M that
% contains cluster number for each observation in the original data.
%
% When bottom leaves are cutoff, some information are lost. T
% supplies this lost information. For example, to find out which
% observations are contained in node k of the dendrogram, use
% FIND(T==k).
%
% When there are less than P observations in the original data, T
% is the identical map, i.e. T = (1:M)'. Each node only contains
% itself.
%
% See also LINKAGE, PDIST, CLUSTER, CLUSTERDATA, INCONSISTENT.
% ZP You, 3-10-98
% Copyright (c) 1993-98 by The MathWorks, Inc.
% $Revision: 1.2 $
m = size(Z,1)+1;
%if nargin < 2
p = 30;
%end
Z = transz(Z); % convert from m+k indexing to min(i,j) indexing.
T = (1:m)';
% if there are more than p node, dendrogram looks crowded, the following code
% will make the last p link nodes as the leaf node.
if (m > p) & (p ~= 0)
Y = Z((m-p+1):end,:);
R = Y(:,1:2);
R = unique(R(:));
Rlp = R(R<=p);
Rgp = R(R>p);
W(Rlp) = Rlp;
W(Rgp) = setdiff(1:p, Rlp);
W = W';
T(R) = W(R);
% computer all the leaf that each node (in the last 30 row) has
for i = 1:p
c = R(i);
T = clusternum(Z,T,W(c),c,m-p+1,0); % assign to it's leaves.
end
Y(:,1) = W(Y(:,1));
Y(:,2) = W(Y(:,2));
Z = Y;
Z(:,3) = Z(:,3)-min(Z(:,3))*0.8; % this is to make the graph look more uniform.
m = p; % reset the number of node to be 30 (row number = 29).
end
A = zeros(4,m-1);
B = A;
n = m;
X = 1:n;
Y = zeros(n,1);
r = Y;
% arrange Z into W so that there will be no crossing in the dendrogram.
W = zeros(size(Z));
W(1,:) = Z(1,:);
nsw = zeros(n,1); rsw = nsw;
nsw(Z(1,1:2)) = 1; rsw(1) = 1;
k = 2; s = 2;
while (k < n)
i = s;
while rsw(i) | ~any(nsw(Z(i,1:2)))
if rsw(i) & i == s, s = s+1; end
i = i+1;
end
W(k,:) = Z(i,:);
nsw(Z(i,1:2)) = 1;
rsw(i) = 1;
if s == i, s = s+1; end
k = k+1;
end
%%%%%% get an X out of W
clear memb c g1 g2
l=length(Z)+1;
%Z=transz(Z1);
c(:,1)=ones(1,l)';
memb(1:l,1)=(1:l)';
a=Z(1,1);b=Z(1,2);
g1(1,1)=a;g2(1,1)=b;
memb(a,1:2)=union(memb(a,1),memb(b,1));
c(a,1)=2;
for k=2:length(Z),
c(:,k)=c(:,k-1);
a=Z(k,1);b=Z(k,2);
g1(k,1:c(a,k))= memb(a,1:c(a,k));
g2(k,1:c(b,k))=memb(b,1:c(b,k));
c1=c(a,k)+c(b,k);
memb(a,1:c1)=union(memb(a,1:c(a,k)),memb(b,1:c(b,k)));
c(a,k)=c1;
end;
Ds=sum(D');
index=1:l;
for k=length(Z):-1:2,
a=Z(k,1);b=Z(k,2);
gr1=g1(k,1:c(a,k-1))
gr2=g2(k,1:c(b,k-1))
t1=mean(Ds(gr1));
t2=mean(Ds(gr2));
if (t1>t2),
[x,y]=sort(index([gr1 gr2]));
index(gr1)=x(1:length(gr1));
index(gr2)=x((length(gr1)+1):end);
end;
if (t1<=t2),
[x,y]=sort(index([gr1 gr2]));
index(gr2)=x(1:length(gr2));
index(gr1)=x((length(gr2)+1):end);
end;
index
end;
%%%%%%%%%%%%%%
X=index;
[u,v]=sort(X);
T=v;
label = num2str(v');
for n=1:(m-1)
i = Z(n,1); j = Z(n,2); w = Z(n,3);
A(:,n) = [X(i) X(i) X(j) X(j)]';
B(:,n) = [Y(i) w w Y(j)]';
X(i) = (X(i)+X(j))/2; Y(i) = w;
end
figure
set(gcf,'Position', [50, 50, 800, 500]);
h = plot(A,B,'b');
axis([0 m+2 0 max(Z(:,3))*1.05])
set(gca,'XTickLabel',[],'XTick',[],'box','off','xcolor','w');
text((1:m)-0.2,zeros(m,1)-0.05*max(Z(:,3)),label);