function [class,CC,prior,postprob,maxcriteria] = CCA(D,Kstart,Kfinal,Rn)
% Usage:
%        [class,CC,prior,postprob,maxcriteria] = CCA(D,Kstart,Kfinal,Rn)
%
% This program computes the clustered components analysis of the data set
% under the assumption that the noise is additive white with convariance
% matrix I. That is, each component is independent with unit variance.
% If this is not the case, the data should be prewhitened to make this so.
%
% input:
%   D(N,M) - data matrix containing N vectors of length M
%   Kstart - initial number of clusters
%   Kfinal - final number of clusters: if Kfinal==0 => estimate # of clusters
%   Rn - Assumed covariance matrix of noise
%        If the noise is i.i.d. with variance SigmaN2, then 
%        Rn = SigmaN2*diag(ones(M,1));
% output:
%   class(N) - hard classification of data vectors
%   CC(M,K) - final values for K cluster component vectors, each of length M
%   prior(K) - prior probabilities for each cluster
%   posterprob(N,K) - posterior probabilities for each data vector and 
%                     each cluster
%   maxcriteria - maximum MDL criteria for final number of clusters
% 
% written by C. A. Bouman
% modified by S. Chen



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This section of code set certain preferred options                             %
%                                                                                %
% This factor can be used to adjust the weighting of the weighting               %
% of the model order term. It is normally set to 1.                              %
%     If OIDfact > 1, this tends to reduce estimated model order.                %
%     If OIDfact < 1, this tends to increase estimated model order.              %
OIDfact = 1;

%   option:                                                                      %
%     option ==0 => (Preferred Option) Treat first cluster the same as others    %
%                   With this option, each vector is assumed                     %
%                   to contain nonzero signal components.                        %
%     option ==1 => Treat first cluster as centered about 0                      %
%                   Use this option if you want to simmultaneously               %
%                   detect signals and estimate their components.                %
option = 0;

%                                                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Compute whitening filter
[U,S,V] = svd(Rn);
lambda = diag(S);
T =  diag(1./sqrt(lambda))*V';
Tinv =  V*diag(sqrt(lambda));

% Whiten data
D = D*T';


% If option==1, treat first class differently
if option~=1
  option = 0;
end

% Seed algorithm 
[CC,prior] = SeedCluster(D,Kstart,option);

% Use EM iterations to estimate CC and prior
[R,CC,prior,postprob,criteria] = EMopt(D,CC,prior,OIDfact,option);

% Compute minimum number of clusters
if option==1
  Kmin = max(Kfinal,2);
else
  Kmin = max(Kfinal,1);
end

% Save state
maxcriteria = criteria;
bestCC = CC;
bestprior = prior;
bestpostprob = postprob;

% Apply iterative cluster merging
for k=Kstart:-1:(Kmin+1)
  % Merge two clusters
  [CC,prior] = MergeCluster(R,prior,option);

  % Use EM iterations to estimate CC and prior
  [R,CC,prior,postprob,criteria] = EMopt(D,CC,prior,OIDfact,option);

  % Save state if it is better then previous
  if criteria > maxcriteria
    maxcriteria = criteria;
    bestCC = CC;
    bestprior = prior;
    bestpostprob = postprob;
  end
end

% Apply order estimation if desired
if Kfinal==0
  criteria = maxcriteria;
  CC = bestCC;
  prior = bestprior;
  postprob = bestpostprob;
  disp(sprintf('Estimated number of clusters = %d',length(prior)));
end

% Reorder classes according to likelihood
tmpprior = prior;
if option==1
  tmpprior(1) = 2;
end
[tmpprior,I] = sort(-tmpprior);
prior = prior(I);
CC = CC(:,I);
postprob = postprob(:,I);

% Compute hard classification
[value,class] = max(postprob,[],2);

% Transform back clustered components vectors
CC = Tinv*CC;
CC = normalize(CC);


