function eta=pemss(z,etainit,index,maxiter,tol,lim,maxsize,T)
%PEMSS   Minimizes prediction error criteria
%
%        TH = pemss(Z,TH0)
%
%	 TH: The result is returned in this matrix. See HELP THETA
%	 Z: The datamatrix Z = [Y U], where each column of Y and U
%	    corresponds to a particular output and input signal.
%	 TH0: This argument provides the model structure,
%	 	the initial conditions, and the quadratic norm.
%		(See MS2TH, MF2TH and THINIT)
%
%	 With TH = pemss(Z,TH0,INDEX) only the adjustable parameters
%	 with indices listed in the row vector INDEX will be estimated
%	 Remaining parameters will be fixed to their values in ETA0
%
%	 Optional arguments can be reached by
%	 TH = pemss(Z,TH0,INDEX,MAXITER,TOL,LIM,MAXSIZE,T)
%	 See Help AUXVAR for an explanation of the these arguments.
%	 In the m-file NUDERSTEP the stepsize to be used for numerical 
%	 differentiation is set.

%        L. Ljung 10-2-90
%        Copyright (c) 1986-92 by the MathWorks, Inc.
%        All Rights Reserved.

[Ncap,nz]=size(z);
maxdef=idmsize(Ncap);
if nargin<8,T=etainit(1,2);end
if nargin<7,maxsize=maxdef;end
if nargin<6,lim=1.7;end
if nargin<5,tol=0.01;end
if nargin<4,maxiter=10;end
if nargin<3,index=[];end
[nr,nc]=size(etainit);
nd=etainit(1,5);
if nd==0,error('This model structure does not have any adjustable parameters!'),end

if isempty(index),index=1:nd;end
sspm=getmfth(etainit);
if any(etainit(2,8)==[2 3 4]), Tmod=-abs(T);else Tmod=abs(T);end
[th,P,lambda]=th2par(etainit);
[modarg,rarg,carg]=getargth(etainit);
n=nd;
[A,B,C,D,K,X0]=feval(sspm,th,Tmod,modarg);[ny,nx]=size(C);
[nx,nu]=size(B);[ny,nx]=size(C);
rowmax=max(length(index)*ny,nx+nz);
Mloop=floor(maxsize/rowmax);
if (nu+ny)~=nz, error('Incorrect number of inputs and outputs specified in data vector'),end
ei=eig(A-K*C);
if max(abs(ei))>1.1,error('Initial guess predictor unstable!'),end
for kc=1:Mloop:Ncap-1
      jj=(kc:min(Ncap,kc-1+Mloop));
      if jj(length(jj))<Ncap,jjz=[jj,jj(length(jj))+1];else jjz=jj;end
      xh=ltitr(A-K*C,[K B-K*D],z(jjz,:),X0);
      yh(jj,:)=(C*xh(1:length(jj),:)'+[zeros(ny,ny) D]*z(jj,:)')';
      X0=xh(length(xh),:)';
end
e=z(:,1:ny)-yh;
V=det(e'*e/length(e));

g=ones(n,1); 
%** Determine limit for robustification of criterion ***

if lim~=0, lim=median(abs(e-ones(Ncap,1)*median(e)))*lim/0.7;end
if lim==0, lim=1000*ones(1,ny)*V;end
        
el=min(e,ones(Ncap,1)*lim); el=max(el,-ones(Ncap,1)*lim);
V=det(el'*e/length(e));
clear e, l=0; st=0;

%  ** The minimization loop **

while [norm(g)>tol l<maxiter st~=1]
   l=l+1;

%   * Compute the Gauss-Newton direction *
   if norm(lambda-eye(ny))==0,lam=eye(ny);else lam=inv(sqrtm(lambda));end
   [g,R,ng]=gnns(z,th,sspm,Tmod,modarg,lam,lim,maxsize,index);

%    * Search along the g-direction for a lower value of V *

    [Vn,th1,st,lambda]=sess(V,z,th,sspm,Tmod,modarg,lam,g,lim,Mloop);

    if st==1,
[Vn,th1,st,lambda]=...
sess(V,z,th,sspm,Tmod,modarg,lam,ng/trace(R)*length(R),lim,Mloop);
    end

%  * Display the current values along with the previous ones *

     t=zeros(n,3); t(:,1)=th1';t(:,2)=th';t(:,3)=g;
     home,clc
     disp(['    ITERATION # ' int2str(l)])
     disp(['Current loss: ' num2str(Vn) '   Previous loss: ' num2str(V)])
      disp(['Current th  prev. th   gn-dir. '])
      disp(t)
      disp(['Norm of gn-vector: ' num2str(norm(g))])
      if st==1
          disp('No improvement of the criterion possible along the search direction')
          disp('Iterations therefore terminated')
      end
      th=th1;V=Vn;
end

% *** Build up the ETA-matrix ***

eta = etainit;
P1=zeros(nd,nd);
P1(index,index)=inv(R);
P=P1;
eta(3,1:nd)=th;
eta(1,1)=det(lambda);eta(2,1)=eta(1,1)*(1+nd/Ncap)/(1-nd/Ncap);
eta(4:nd+3,1:nd)=P;
eta(4+nd+etainit(1,6):3+ny+nd+etainit(1,6),1:ny)=lambda;
ti=fix(clock); ti(1)=ti(1)/100;
eta(2,2:6)=ti(1:5); eta(2,7)=23;
if eta(2,8)==3,eta(2,7)=33;end

