function Clusters=hierclst(vect,k)
% HIERCLST  Hierarchical clustering
%	Clusters=hierclst(vect,k) return the cluster pairs
%	when the dataset vect (matrix of column vectors) is clustered
%	first into k clusters, then each cluster is clustered into k clusters,
%	etc. The returned matrix has 1's indicating cluster 1, 2's, 3's etc.
%	k-means clustering is used at each step.
%	
%	If k is a vector, the k-means clustering at level i uses k=k(i).
%	If the number of levels exceeds the length of k, the last element
%	of k is used for all the lower levels.

k=[k k(length(k))];
if size(vect,2)>k(1)
  [dummy,Clusters]=kmeans(vect,k(1));
  AllNewClusters=[];
  for i=1:k(1)
    Idx=find(Clusters==i);
    TheseNewClusters=hierclst(vect(:,Idx),k(2:length(k)));
    if size(TheseNewClusters,1)>size(AllNewClusters,1)
      % Resize AllNewClusters
      AllNewClusters(size(TheseNewClusters,1),size(Clusters,2))=0;
    end
    AllNewClusters(1:size(TheseNewClusters,1),Idx)=TheseNewClusters;
  end
  Clusters=[Clusters;AllNewClusters];
else
  Clusters=1:size(vect,2);
end