function th=oe(z,nn,maxiter,tol,lim,maxsize,Tsamp)
%OE	Computes the prediction error estimate of an output-error model.
%
%	TH=oe(Z,NN)
%
%	TH: Returned as the estimated parameters of the output-error model
%	y(t) = [B(q)/F(q)] u(t-nk) + e(t)
%	along with estimated covariances and structure information.
%	For the exact format of TH, see HELP THETA.
%
%	Z : The output-input data Z=[y u], with y and u being column vectors.
%	This routine only works for single-input systems. Use PEM otherwise.
%
%	NN: Initial value and structure information, given either as
%	NN=[nb nf nk], the orders and delay of the above model, or as
%	NN=THI, with THI being a theta matrix of the standrad format. In
%	the latter case the criterion minimization is initialized at THI.
%
%	Some parameters associated with the algorithm are accessed by
%	TH = oe(Z,NN,maxiter,tol,lim,maxsize,T)
%	See HELP AUXVAR for an explanation of these, and their default values.

%	L. Ljung 10-1-86,4-7-92
%	Copyright (c) 1986-92 by the MathWorks, Inc.
%	All Rights Reserved.

% *** Set up default values ***
[nr,cc]=size(nn); [Ncap,nz]=size(z); nu=nz-1;
maxsdef=idmsize(Ncap);

if nargin<7, Tsamp=[];end
if nargin<6, maxsize=[];end
if nargin<5, lim=[];end
if nargin<4, tol=[];, end
if nargin<3, maxiter=[]; end
if isempty(Tsamp),if nr>1,Tsamp=gett(nn);else Tsamp=1;end,end
if isempty(maxsize),maxsize=maxsdef;end
if isempty(lim),lim=1.6;end,if isempty(tol),tol=0.01;end
if isempty(maxiter),maxiter=10;end
if Tsamp<0,Tsamp=1;end, if maxsize<0,maxsize=maxsdef;end, if lim<0,lim=1.6;end
if tol<0,tol=0.01;end,if maxiter<0, maxiter=10;end


%  *** Do some consistency tests ***

if nz>Ncap, error('The data should be organized in column vectors')
return,end
if nu>1, error('This routine only works for single-input systems. Use PEM instead!')
return,end
if nu==0, error('OE makes sense only if an input is present!')
return,end
if nr==1 & cc~=3,
         error('Incorrect number of orders specified: Should be  nn=[nb nf nk]')
return,end
%
%*** if nn has more than 1 row, then it is a theta matrix
%     and we jump directly to the minimization ***
%
if nr>1 nu=nn(1,3);
   if nu>1,
      error('This routine only works for single-input systems')
   return,end
   nf=nn(1,8); nb=nn(1,5); nk=nn(1,9); n=nf+nb;ni=max(nf,nb+nk-1);
   t=nn(3,1:n)';
   if nf>0; f=[1 t(nb+1:n)'];else f=1;end
end
%
%
if nr==1
   nb=nn(1); nf=nn(2); nk=nn(3); n=nb+nf;ni=max(nf,nb+nk-1);
%
%
% *** compute initial estimate via LS and IV ***
%
  t=arx(z,[nf nb nk],maxsize,Tsamp,0);
  f=fstab([1 t(1:nf)']);
        
  t=iv(z,[nf nb nk],f,[zeros(1,nk) t(nf+1:n)'],maxsize,Tsamp,0);
  f=fstab([1 t(1:nf)']);
        
%
  if nf>0, t=[t(nf+1:n); t(1:nf)];end
end
% *** Display initial estimate ***
%
w=pefilt([zeros(1,nk) t(1:nb)'],f,z(:,2),z(1:ni,1));
v=z(:,1)-w; Vcap=v'*v/(Ncap-ni);
if isnan(Vcap)|~finite(Vcap)
	error('Your data or model structure has some deficiencies. Try a different input or a lower order model!'),end
clc, disp(['  INITIAL ESTIMATE'])
disp(['Current loss:' num2str(Vcap)])
disp(['th-vector'])
disp(t)
%
% *** start minimizing ***
%
% ** determine limit for robustification **
 if lim~=0, lim=median(abs(v-median(v)))*lim/0.7;end
 if lim==0, lim=1000*Vcap;end
 [mmv,nnv]=size(v);
 ll=ones(mmv,nnv)*lim;vl=max(min(v,ll),-ll); clear ll
 g=ones(n,1); l=0; st=0; nmax=max(nf,nb+nk-1);
%
% ** the minimization loop **
%
while [norm(g)>tol l<maxiter st~=1]
      l=l+1;
%     * compute gradient *
      wf=filter(-1,f,w); uf=filter(1,f,z(:,2));
   M=floor(maxsize/n);
   R=zeros(n);Fcap=zeros(n,1);
   for k=nmax:M:Ncap-1
      jj=(k+1:min(Ncap,k+M));
      psi=zeros(length(jj),n);
      for k=1:nb, psi(:,k)=uf(jj-k-nk+1);end
      for k=1:nf, psi(:,k+nb)=wf(jj-k);end
      if Ncap>M,R=R+psi'*psi; Fcap=Fcap+psi'*vl(jj);end
   end
   if Ncap>M, g=R\Fcap; else g=psi\vl(jj);end
%
%     * search along the g-direction *
%
      [t1,w,vl,v,V1,f,st]=searchoe(z,t,g,lim,Vcap,nb,nf,nk,ni);
      home
      disp(['      ITERATION # ' int2str(l)])
      disp(['Current loss: ' num2str(V1) '  Previous loss: ' num2str(Vcap)])
      disp(['Current th  prev. th   gn-dir'])
      disp([t1 t g])
      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
      t=t1; Vcap=V1;
end
Vcap=v'*v/(length(v)-ni);
th(1,:)=[Vcap Tsamp 1 0 nb 0 0 nf nk];
th(2,1)=Vcap*(1+n/Ncap)/(1-n/Ncap);
ti=fix(clock); ti(1)=ti(1)/100;
th(2,2:6)=ti(1:5);
th(2,7)=8;
th(3,1:n)=t';
if maxiter==0, return,end
if Ncap>M, PP=inv(R); else PP=inv(psi'*psi);end
th(4:3+n,1:n)=Vcap*PP;

