function [par,struc,it_inf,cov,lambda] = minloop(z,par,struc,algorithm,it_inf)
%MINLOOP  ** The minimization loop **

%   Copyright 1986-2000 The MathWorks, Inc.
%   $Revision: 1.14 $  $Date: 2000/09/01 13:41:19 $

Xsum = findobj(get(0,'children'),'tag','sitb16');
if ~isempty(Xsum)
   XID = get(Xsum,'UserData');
   try
      XIDiter = XID.iter;
   catch
   end
end
oeflag = 0;
switch struc.type
case 'poly'
   if struc.na+struc.nc+struc.nd==0
      oeflag = 1;
   end
case 'ssnans'
   if ~any(any(isnan(struc.filearg.ks))')
      if norm(struc.filearg.ks)==0
         oeflag = 1;
      end
   end
case 'ssgen'
   m0 =struc.model;
   lnp = length(pvget(m0,'ParameterVector'));
   m0 = parset(m0,randn(lnp,1));
   if norm(pvget(m0,'K'))==0
      oeflag = 1;
   end
   
end
algorithm.FixedParameter = pnam2num(algorithm.FixedParameter,struc.pname);
if (length(algorithm.FixedParameter) == length(par))&length(par)>0
   error('All parameters are fixed. Nothing to estimate.')
end

gui=0;
stopp=0;flag=[];
try
flag=findobj(get(0,'children'),'flat','tag','sitb22');
end
if ~isempty(flag),
   if strcmp(get(flag,'vis'),'on'),
      gui=1;set(XIDiter(7),'userdata',0);
   end,
end 
npar = length(par); % FUnkar ej vid free!!!
tol = algorithm.Tolerance;
maxiter = abs(algorithm.MaxIter);
impr=[];
ind1 = algorithm.FixedParameter;
ind2 = [1:npar]';
if isempty(ind1)
   ind3=ind2;
else
   indt = ind2*ones(1,length(ind1))~=ones(npar,1)*ind1(:)';
   if size(indt,2)>1, indt=all(indt');end
   ind3 = ind2(find(indt));%all((ind2*ones(1,length(ind1))~=ones(npar,1)*ind1)')))
end
%%LL FIXA LIM
algorithm.estindex=ind3;
lambda=struc.lambda;
if norm(lambda==0) %% This is to protect from strange initial model
   lambda=eye(size(lambda));
end

ny=1;                  
testnorm = tol+1;
stop = 0;stopp=0;
l = 0;
if strcmp(lower(algorithm.Trace),'on')
   dispmode = 1;
elseif strcmp(lower(algorithm.Trace),'full')
   dispmode = 2;
else
   dispmode = 0;
end

while [testnorm>tol l<maxiter stop~=1]
   l=l+1;
   
   if gui,
      drawnow,
      display=get(XIDiter(2),'value');
      if display,dispmode = 2;end
   end
   % Compute the Gauss-Newton direction *
   struc.sqrlam=inv(sqrtm(lambda));
   
   if strcmp(struc.type,'ssfree')
      if ~isempty(par)
         struc=ssfrupd(par,struc);
      end
      par=[];
      dkx=struc.dkx;m=struc.nu;p=struc.ny;n=struc.nx;nk=struc.nk;
      if dkx(2), bk = [struc.b struc.k];, else, bk = struc.b; end
      Qperp = fipert_qr(struc.a,bk,struc.c);
      if isempty(nk)
         nextra = dkx(3)*n;
      else
         nextra = dkx(1)*sum(nk==0)*p+dkx(3)*n;
      end
      if nextra>0,
         Qperp = [Qperp zeros(size(Qperp,1),nextra); zeros(nextra, ...
               size(Qperp,2)) eye(nextra)];
      end
      struc.Qperp=Qperp;
      
   end
   [lambda,tlam,psi,e,Nobs]=gnnew(z,par,struc,algorithm);V=det(lambda);    
   if isempty(psi)
      error('Predictor has become unstable')
   end
   if algorithm.MaxIter<0
      parnew=par;
      lambda = tlam;
      break
   end
   
   if l==1
      delta = algorithm.Advanced.Search.LmStartValue*norm(psi);
   end
   testnorm=abs((e'*psi)*pinv(psi'*psi)*(psi'*e)*100/Nobs);  
   % No normalization wrt V for ss, since that   
   % is implicit in the lam-normalization 
   if strcmp(struc.type,'poly')
      testnorm=testnorm/lambda;
   end
   if dispmode
      disp(['Iteration # ',int2str(l),':'])
   end
   
   [Vn,parnew,stop,delta,lambda,tlam]=msearch(V,z,par,struc,psi,e,delta,algorithm,...
      lambda,tlam,dispmode); 
   %  * Display the current values along with the previous ones *
   
   g=pinv(psi)*e;
   
   if dispmode      
      disp(['   Current loss: ' num2str(Vn) '   Previous loss: ' num2str(V)])
      if ~strcmp(struc.type,'ssfree')&dispmode>1
         t=zeros(npar,3); t(:,1)=parnew;t(:,2)=par;
         t(algorithm.estindex,3)=pinv(psi)*e;
         disp(['   New par   prev. par   gn-dir '])
         disp(t)
      end
      disp(['   Norm of gn-vector: ' num2str(norm(g))])
      disp(['   Expected improvement: ',num2str(testnorm),' %'])
      disp(['   Achieved improvement: ',num2str((V-Vn)*100/V),' %'])
      if stop==1
         disp('   No improvement of the criterion possible along the search direction.')
         disp('   Iterations therefore terminated.')
      end
   end %if display
   if gui
      set(XIDiter(4),'string',int2str(l));
      set(XIDiter(5),'string',num2str(Vn));
      set(XIDiter(9),'string',['(',num2str(V),')']);
      set(XIDiter(6),'string',num2str(testnorm));
      drawnow
      stopp=get(XIDiter(7),'userdata');
   end
   testnorm = max(testnorm,(V-Vn)*100/V); 
   par = parnew;impr = -(Vn-V)*100/V;V=Vn;
   V = det(tlam); 
   lambda = tlam;% For the noise Cov estimation we do not use robustification
   if strcmp(struc.type,'ssfree')
      struc=ssfrupd(par,struc);
      par=[];
   end
   if stopp,break,end
end
if l==maxiter
   it_inf.WhyStop='Maxiter reached';
elseif stop
   it_inf.WhyStop='No improvement along search direction.';
else
   it_inf.WhyStop='Near (local) minimum, (norm(g)<tol).';
end
if strcmp(struc.type,'ssfree')
   npar = length(psi);
end

it_inf.UpdateNorm=testnorm;
it_inf.LastImprovement=[num2str(impr),'%']; 
it_inf.Iterations=l;
Vl = V*(1-npar/Nobs); % V has been compensated for # of pars
it_inf.LossFcn=Vl;
it_inf.FPE = Vl*(1+npar/Nobs)/(1-npar/Nobs);
it_inf.InitialState = struc.init;
if strcmp(struc.type,'ssfree')
   cov=[];
elseif ~struc.cov
   cov = 'None';
else
   cov = zeros(npar,npar);
   pin=pinv(psi'*psi); 
   if oeflag    
      [lambda,tlam,psi,e,Nobs]=gnnew(z,par,struc,algorithm,1);
      pin = pin*(psi'*psi)*pin;
      lambda = tlam; % For the noise Cov estimation we do not use robustification
   end
   cov(ind3,ind3) = pin;
   if strcmp(struc.type,'poly')&~oeflag % the psi is not normalized by lamdba
      cov = V*cov;                
   end
   
   try
      norm(cov);
   catch
      disp('Covariance Matrix Ill-conditioned. Not stored.')
      cov = [];
   end
end


