computeProportionalModel(covList, m_constants)

function [proportionalCovList, LList, B, rhoVec] = computeProportionalModel(covList, m_constants)

% Proportional Covariance Matrices Algorithm

% Flury, Bernhard. Commom Principal Components and Related multivariate
% Models, Wiley
    
    
     [cpcCovList, LList, B, FList] = computeCPCModel(covList, m_constants);
     if isempty(cpcCovList)
         rhoVec = ones(m_constants.numGrp,1)'; % Usual initial approximation
     else
        rhoVec(1) = 1;
	    V1 = cpcCovList{1};
	    for i = 2:m_constants.numGrp 
            rhoVec(i) = (det(cpcCovList{i}) / det(V1))^(1 / m_constants.numGrp);
	    end  
    end

     rhoVecOld = rhoVec;
     epsf = 10.^(-5);
     diff = 1;
     coefMLElambda = (m_constants.aprioriProb ./ rhoVec);
     while diff > 0
         tempCov = zeros(m_constants.numVar);
         for i = 1:m_constants.numGrp 
            tempCov = tempCov + coefMLElambda(i) * covList{i};
         end
         [B,L] = eig(tempCov);
         lambdaMat = diag(B' * covList{1}' * B)';
         for i = 2:m_constants.numGrp 
            lambdaMat = [lambdaMat;diag(B' * covList{i}' * B)'];
         end
         lambdaVec = coefMLElambda * lambdaMat;
         coefMLErho = ((m_constants.numVar * lambdaVec).^-1);
         rhoVec(2:end) = lambdaMat(2:end,:) * coefMLErho';
         diffRho = abs(rhoVec - rhoVecOld);
         eps = max(diffRho);
         diff = eps - epsf;
         rhoVecOld = rhoVec;
         coefMLElambda = (m_constants.aprioriProb ./ rhoVec);
     end
     for i = 1:m_constants.numGrp
        LList{i} = rhoVec(i) * lambdaVec;
        proportionalCovList{i} = B * diag(LList{i}) * B';
	end