root/library/utia_legacy/ticket_37/uoptim.m @ 653

function [aMix] = soptim(aMix, aMixu, ufc, nstep, chis)
% soptim performs simultaneous advisory design for normal mixture
%
% [aMix] = soptim(aMix, aMixu, ufc, nstep, chis) 
% [aMix] = soptim(aMix, aMixu, ufc, nstep) chis = 1 
% [aMix] = soptim(aMix, aMixu, ufc)        nstep = [200, 1] 
% 
% aMix   : advised mixture of the type ARX LS enriched on following control states:
%            strc : common control structure 
%            ufc  : normalised vector qualifying components:  
%                   dangerous component (0), not dangerous (positive number) 
%            kc   : lift of quadratic forms            
%            UDc  : cell vector of u'du decompositions of KLD kernels 
%            udca : u'du decomposition of average KLD kernel in UDc 
%            kca  : average lift of quadratic forms kc
% aMixu  : desired mixture (user's target) of the type ARX LS with control states
% ufc    : vector qualifying components: 0 - dangerous component, (1) - not 
% nstep  : parameters [ns1,per] determining design horizon, i.e. horizon = ns1*per; 
%            ns1  : number of block repetition
%            per  : horizon of a block  
%          if nstep is defined by parameter nsl only then per is set to 1 
% chis   : indicates strategy chosen: chis=1 for receding horizon (default) and chis=-1 for IST
%
% Design     : J. Bohm
% Updated    : June, 2002
% Project    : ProDaCTools, IST-1999-12058 
% See also   : udupdt, getdvect, facchng, facarxls

% References : \ref{ch9} 
% Note       : 
% Updated    :  


if (nargin < 3) | ~any(ufc)error(sprintf('%s\n%s\n%s','ufc was not correcly set, define it as a vector','of a length "ncom" having at least one nonzero element ','or use a function ufcgen')) ;end

if (nargin <4), nstep=[200,1];chis=1;end
if (nargin <5), chis=1;end

if(isempty(aMix.states.uchn)) error('uoptim needs nonempty list of channels with recognisable actions'); end;
% normalisation of ufc
ufc = ufc/sum(ufc);

%Inititialization
  ncom   = length(aMix.dfcs);   
  dfcs   = aMix.dfcs;
  strc   = aMix.states.strc;             % common control structure
  nPsi   = max(size(strc));              % length of regression vector + data 
  pochn  = aMix.states.pochn;            % list of channels with o-innovations
  npochn = length(pochn);                % number of channels with o-innovations
  nychn  = length(aMix.states.modelled); % number of modelled channels
  nouts  = length(aMix.states.outs);     % number of innovation channels
  npsi   = nPsi-nychn;                   % length of the regression vector  
  kc0    = aMixu.states.kc;               
  udca   = aMixu.states.udca;
  kca    = aMixu.states.kca;
  coms   = aMix.coms;
  Ethz   = zeros(1,nPsi);
  UDc    = aMixu.states.UDc;
  lss    = length(nstep);
       lrica  = zeros(1,nPsi-1);
  %test of aMixu  
 coves=zeros(1,npochn);
 % if npochn~=length(aMixu.Facs),error('aMixu not correctly set'); end
  for i=1:npochn,
          coves(i)=aMixu.Facs{i}.cove;
  end
  if ~any(coves), error('aMixu not correctly set, Facs{.}.cove must be >0'),end

   for i=1:ncom,                         % cycle over number of components ncom
          Ric{i}= zeros(nPsi);               % KLD kernels 
          kcc(i)=0;                          % lift of quadratic forms
   end  
   df=dfcs/sum(dfcs);
 % setting of design horizon  
  if lss==2 
      steps  = nstep(1)*nstep(2); per=nstep(2);
  else 
      steps = nstep; per=1;
  end  
  if chis>0,                            % if the strategy starts from zero
         udca = zeros(nPsi); kca =0;
  end  
 %Main design cycle, iterations over the horizon of the criterion 
  for iter=1:steps, % =====================  iterations till design horizon
%  ricmn is an auxiliary array accumulating  results of optimization  
      if mod(iter-1,per)==0,
% shift of a matrix from bottom right to top left by nychn  
  %       if nPsi>nychn+1
            [udca, lrica]= ricshift(udca,lrica,nychn,nPsi,npsi);   
  %     end %if nPsi
         udca(nPsi,nPsi)= 0;    
     end %if mod
 
     for i=1:ncom % ......................... cycle over all components 
         if mod(iter-1,per)==0,  
   %         if ufc(i)==0,  kcc(i)=1e30; continue; end           % excluding bad components
                ric = udca;    % ric is auxiliarry working array
            lric(i,:)=lrica;
            kcc(i) = -npochn +kc0(i);  
 % adding to each component its stationary loss       
           for j=1:npsi,   % 
               red=UDc{i}(j,:);
               red(j)=1;
               ric= udupdt(ric,red,ufc(i)*UDc{i}(j,j));
           end  % for j
        else  
% iteration continues in corresponding component kernel and lift
% shift of a matrix from bottom right to top left by nychn 
           if nPsi>nychn+1, 
               [ric,lric(i,:)]=ricshift(Ric{i},lric(i,:),nychn,nPsi,npsi);
           end % if nPsi
           ric(nPsi,nPsi) = Ric{i}(nPsi,nPsi);
        end   %if  mod   
% expectation is calculated channel by channel   
          for j=1: nouts %---------------------  cycle over innovation  channels
                      indv = ~isempty(find(strc(1,j)==pochn));    % indicator if the channel is o-innovation
              [ric,lric,kcc]= ricexp(ric,lric,kcc,i,j,aMix,nPsi); 
        % the penalization is used 
              if indv,                           % visibility indicator
             [ric,lric,kcc]= ricpen(ric,lric,kcc,i,j,aMix,aMixu,nPsi);  
              end   %end if indv
         end % ------------------------------ reduced all factors of i-th component
  % now penalization
              for j=nouts+1: nychn,
                 [ric,lric,kcc]= ricpenu(ric,lric,kcc,i,j,aMix,aMixu,nPsi); 
  %               if ric(j,j)>eps, 
  %                  kcc(i)=kcc(i)-lric(i,j)*lric(i,j)/ric(j,j)/4;
  %                  disp('pred m');
  %                  keyboard
% %                   lric(i,:)=lric(i,:)-lric(i,j)*[zeros(1,j) ric(j,j+1:end-1)];
  %                  disp('po  m');
  %                  keyboard
  %             end  % if  
         end   %end for
       
        [l,d]= ld2ld(ric(nouts+1: nychn,nouts+1: nychn));
           lric(i,nychn+1:end)=lric(i,nychn+1:end)-lric(i,nouts+1: nychn)*inv(l)*ric(nouts+1: nychn,nychn+1:end-1);
             Ric{i} = ric;
     end          % ...................... done for all components
% now putting losses together
         if mod(iter,per)==0,   
            udca=zeros(size(ric));
           lrica  = zeros(1,nPsi-1);
            det=ones(1,ncom);
%%putting it together   
            for i=1:ncom,
               % if ufc(i)==0, continue; end
                for j=nouts+1:nPsi,   
                        red = Ric{i}(j,:);
                        red(j) = 1;
                    udca = udupdt(udca,red,df(i)*Ric{i}(j,j));
             end       % over j 
              lrica=lrica+df(i)*lric(i,:);

            end          % over components
    end   %if  
      
end            % iterations       
%recalculating lric into the triangular matrix
%   for i=1:ncom % ......................... cycle over all components  
 %     ric=Ric{i};
      [r,d]=ld2ld(udca);
      pom=lrica';
      xx=zeros(nPsi -1,1);
      if d(nouts+1,nouts+1)>eps, xx(nouts+1)=pom(nouts+1)/d(nouts+1,nouts+1)/2;
      else xx(nouts+1)=0; 
      end
      for j=nouts+2:nPsi-1,
          ff=r(nouts+1:end-1,j)'*d(nouts+1:end-1,nouts+1:end-1)*2*xx(nouts+1:end);
          if d(j,j)>eps,xx(j)=(pom(j)-ff)/2/d(j,j);
          else xx(j)=0;
          end
      end 
     udca(:,end)= [xx ;0];  
     %  end     

%puting calculated control factors into  a mixture
       for i=1:ncom,
  %     cl= Ric{i}(nouts+1:nychn,:);
       for j= nouts+1:nychn,  
           jj=aMix.states.strc(1,j);
           Fac=facarxls(jj,strc(:,j+1:end));
           Fac.Eth=-udca(j,j+1:end);
           if udca(j,j)<eps, error('synthesis uoptim uncorrect, check its input parameters'),end
           Fac.cove=1/udca(j,j);
           aMix=facchng(aMix,i,Fac);
       end
   end
     aMix.states.udca = udca;
     aMix.states.kca = kca;
     aMix.states.UDc = Ric;
     aMix.states.kc = kcc;