function [PCC,PCI,est] = pcascal(O,rp,r,coefom,coefa,visualize,corre);
%
% function returns r principal components of matrix O
%
%  rp          - number of components used for estimation of mean value model
%  coefom      - ratio of weights for different models of noise covariance
%                            equal noise model
%                 coefom= ----------------------
%                         independent noise model
%
%                 0 = noise covariance is detected independently
%                 <inf = noise covariance is identity matrix
%
%  coefa       - stabilizing coeficient for estimation of model parametres
%                default (0.01)
%  corre			- use correlations 0=no, 1=yes
global ui_action

%%%%%%%%%%%%%%%%
% user defined 
if nargin<6
  visualize=0;
  if nargin<5
    coefa=0.01;
    if nargin<4
       coefom=1;
       if nargin<3
          rp=r;
       end
    end
 end
end

%%%%%%%%%%%%%%%%
% initialization

T	=size(O,2);
N	=size(O,1);
It	=eye(T);
Ir	=diag([ones(r,1); zeros(T-r,1)]);
Irp	=diag([ones(rp,1); zeros(T-rp,1)]);

%%%%%%%%%%%%%%%%%
% prior analysis

OO	=O'*O;

OO	=OO/trace(OO)*T;%+ones(size(OO));

est	=zeros(3,1); 		% alpha, beta, eps, eva1/eva2, a, b,

beta	=0;
if corre
   Omega=diag(1./diag(OO));
else
   Omega	=It;
end
a	=zeros(T-1,1);
eps	=1;
Delta	=It+diag(a,-1);
Delta(T,T)=eps;
DD	=Delta*Delta';

i=1;

while (i>0 ) & (abs(est(1,i)-beta)>=0.01*beta)
if i==1
  r=2;
else
  r=3;
end;
Ir	=diag([ones(r,1); zeros(T-r,1)]);

   if 1
      C	=Omega+beta*DD;   % pouziti obou modelu SUM + DATA
   else
      C	=Omega;				% pouziti pouze modelu sumu - skalovani
   end      
   H	=chol(C)';
   iH	=It/H;
   FF	=Omega*iH';

   %%%%%%%%%%%%%%%%%%%%%%%%%
   %% principal components
   [EVE,EVA]=eig(FF'*OO*FF);
   
   % re-order eigenvalues and eigenvectors
   EVA=diag(EVA);
   [EVA,e]=sort(EVA);
   EVE=EVE(:,e);
   %
   EVA=flipud(EVA);
   EVE=fliplr(EVE);
   %
   dEVA=EVA;
   EVA=diag(dEVA);
   
   % matrices for estimates 
   %
   pom	=iH'*EVE*(diag(dEVA)*Ir)*EVE';
   M	=pom*iH;
   pomp	=iH'*EVE*(diag(dEVA)*Irp)*EVE';
   Mp	=pomp*iH;
   Q	=pom*H'/Omega;
   P	=M+OO-2*Q;
   dP	=diag(P);

   est(1,i+1)	=beta;
   est(2,i+1)	=mean(a);
   est(3,i+1)	=EVA(1,1)/EVA(2,2);
   if i>1
      e1_e2=est(3,i+1)/est(3,2);
      if (e1_e2>1.1) | (e1_e2)<0.9
         disp('******** Something wrong !! ');
         procent=floor(100*e1_e2);
         disp(['   ratio of eigenvalues  differs by ' num2str(procent) '% !']);
      end;
   end;
   
   %%%%%%%%%%%%%%%%
   %%%% estimates
   mean_dP=mean(dP);
   om=(N+coefom*N)*ones(T,1)./(dP+coefom*mean_dP);
%   om=(N)*ones(T,1)./(dP+coefom*mean_dP);
   Omega=diag(om);

   if 1  %odhad 
     dM=diag(M);
     dM1=diag(M,1);
     
     stabil=mean(dM)*coefa;
     a=-(dM1+stabil)./(dM(1:T-1)+stabil);
     b=-sum(dM1)/sum(dM(2:T))  * ones(T-1,1);
     Delta=eye(T)+diag(a,1); Delta(1,1)=0;
  end

  
   trDDM=trace(Delta*Delta'*Mp);
   beta	= (T-1)*r/(trDDM);%+(T-1)*r/N/20*mean_dP)
   eps	=sqrt(r/beta/(M(1,1)));
        
%   Delta(T,T)=eps;  % dodatecne eps
    Delta(1,1)=eps;  % dodatecne eps
   DD	=Delta*Delta';
   
if i==1
  ipmi_save
end;

    
   if mod(i,visualize)==0
      if i/visualize<1.9
        F_size=get(0,'ScreenSize');
        figure('Position',F_size);
        ui_but_pos=(F_size(4))-80;
        uicontrol('Style','pushbutton','Position',[10,ui_but_pos,60,20],...
                  'String','Finish','Callback',...
                  'global ui_action,ui_action=1;clear ui_action;uiresume');
        uicontrol('Style','pushbutton','Position',[10,ui_but_pos-30,60,20],...
                  'String','Keyboard','Callback',...
                  'global ui_action,ui_action=2;clear ui_action;uiresume');
        uicontrol('Style','pushbutton','Position',[10,ui_but_pos-60,60,20],...
                  'String','Continue','Callback',...
                  'global ui_action,ui_action=3;clear ui_action;uiresume');
      end
      subplot(2,4,1);
      if i>2
        hold off;
        s_e =est(:,3:i+1)./(est(:,3)*ones(1,i-1));
        plot(3:i+1,s_e);
        text(i,s_e(1,i-1),'\beta'); 
        text(i,s_e(2,i-1),'\epsilon'); 
        text(i,s_e(3,i-1),'e1/e2'); 
      end  
      if 1
      subplot(2,4,2); plot(1:T,N*1./om,'r'); title('\Omega');
      subplot(2,4,3); surfc(iH); title('iH');
      subplot(2,4,4); dEVA=abs(diag(EVA)); bar(dEVA(1:3*r));

      subplot(2,4,5); surfc(OO); title('OO');  	ax_OO=axis;
      subplot(2,4,6); surfc(Mp); title('M'); 		axis(ax_OO); % stejne skalovani
      subplot(2,4,7); surfc(r/beta*It/(DD)); 
                      title('\mu_0 = \beta{}DD^{-1}'); %axis(ax_OO); % stejne skalovani
      subplot(2,4,8); plot(1:T-1,a,1:T-1,b);
      end;
      uiwait;
      global ui_action
      if ui_action==1
        break;
      elseif ui_action==2
        keyboard;
      end;
      clear uiaction
   end;

   i=i+1;
   disp(['i=' num2str(i)]);   
end;



ipmi_save

Cur=sqrt(EVA*Irp)*EVE'*iH;
%keyboard
%Cur=EVE'*iH;
PCC=Cur(1:rp,:)';
% pokus o srovnani na kladne hodnoty
PCC=PCC*diag(sign(mean(PCC)));

fid=fopen('matrix1','wb');
fwrite(fid,FF,'double');
fclose(fid);
fid=fopen('matrix2','wb');
fwrite(fid,iH,'double');
fclose(fid);
save matrix1.txt FF -ascii
save matrix2.txt iH -ascii

PCI=O*Omega*PCC/(PCC'*Omega*PCC);
%PCI=O*PCC/(PCC'*PCC);