function V=ivstruc(z,zv,nn,p,maxsize)
%IVSTRUC	Computes the output error fit for families of single-output
%	ARX-models.
%
%	V = ivstruc(ZE,ZV,NN)
%
%	V: The first row of V is returned as the output error fit for 
%	models of the structures, defined by the rows of NN.
%	These models are computed by the IV-method. The next rows
%	of V are returned as NN'. The last row of V contains the logarithms 
%	of the conditioning numbers of the IV-matrix for the structures in 
%	question. The last column of V contains the number of data points in Z
%
%	ZE : the output-input data Z=[y u], with y and u as column vectors.This
%	data is used for estimation. For multi-input systems u=[u1 u2 ... un]
%	ZV : output-input data used to evaluate the models (could equal ZE).
%	NN: Each row of NN is of the format [na nb nk], the orders and 
%	delays associated with the corresponding model. (See HELP IV4 and 
%	HELP STRUC.)
%
%	See HELP SELSTRUC for analysis of V. With
%	V = ivstruc(Z,NN,p,maxsize) the calculation of conditioning numbers
%	can be suppressed (p=0), and the speed/memory trade-off affected.
%	See HELP AUXVAR.

%	L. Ljung 4-12-87,12-27-91
%	Copyright (c) 1987-92 by the MathWorks, Inc.
%	All Rights Reserved.

[Ncap,nz]=size(z); nu=nz-1;[nm,nl]=size(nn);
na=nn(:,1);nb=nn(:,2:1+nu);nk=nn(:,2+nu:1+2*nu);
[snbr,snbc]=size(nb);if any(any(nb==zeros(snbr,snbc))'),
error('nb=0 is not feasible for ivstruc!'),end % Mod 911227
nma=max(na);nbkm=max(nb+nk)-ones(1,nu);nkm=min(nk);
n=nma+sum((nbkm-nkm))+nu;nbm=max(nb);
% 
% *** Some initial tests on the input arguments ***
%
maxsdef=idmsize(Ncap,n);

if nargin<5,maxsize=maxsdef;,end
if nargin<4, p=1;,end

% *** construct instruments (see (7.111)-(7.112)) ***
tha=arx(z,[nma nbm nkm],maxsize,1,0);
NF=fstab([1 tha(1:nma)']);
 	
     x=zeros(Ncap,1);rs=nma;
     for k=1:nu
     	 MF=[zeros(1,nkm(k)) tha(rs+1:rs+nbm(k))'];rs=rs+nbm(k);
         x=x+filter(MF,NF,z(:,k+1));
     end
%
% construct regression matrix
%
nmax=max(max([na+ones(nm,1) nb+nk]'))-1;
M=floor(maxsize/n);nnuu=sum(nbkm)-sum(nkm)+nu;
Rxx=zeros(nma);Ruu=zeros(nnuu);Rxu=zeros(nma,nnuu);Rxy=zeros(nma);
Ruy=zeros(nnuu,nma); F=zeros(n,1);
for k=nmax:M:Ncap
     jj=(k+1:min(Ncap,k+M));
     phix=zeros(length(jj),nma); phiy=phix; phiu=zeros(length(jj),nnuu);
     for kl=1:nma, phiy(:,kl)=-z(jj-kl,1); phix(:,kl)=-x(jj-kl); end
     ss=0;
     for ku=1:nu
            for kl=nkm(ku):nbkm(ku), phiu(:,ss+kl+1-nkm(ku))=z(jj-kl,ku+1);end
             ss=ss+nbkm(ku)-nkm(ku)+1;
     end
     Rxy=Rxy+phix'*phiy;
     Ruy=Ruy+phiu'*phiy;
     Rxu=Rxu+phix'*phiu;
     Ruu=Ruu+phiu'*phiu;
     Rxx=Rxx+phix'*phix;
     if nma>0, F(1:nma)=F(1:nma)+phix'*z(jj,1);end
      F(nma+1:n)=F(nma+1:n)+phiu'*z(jj,1);
end
clear phiu, clear phix, clear phiy, 
%
% compute estimate

for j=1:nm
	s=[];rs=0;
	for ku=1:nu
		s=[s,rs+nk(j,ku)-nkm(ku)+1:rs+nb(j,ku)+nk(j,ku)-nkm(ku)];
		rs=rs+nbkm(ku)-nkm(ku)+1;
	end
	RR=[Rxy(1:na(j),1:na(j)) Rxu(1:na(j),s);Ruy(s,1:na(j)) Ruu(s,s)];
	TH=RR\F([1:na(j) nma+s]);
	a=TH(1:na(j))';b=TH(na(j)+1:length(TH))';f=fstab([1 a]);
	[mzv,nzv] = size(zv(:,1));
	yhat=zeros(mzv,nzv);rs=0;
	for ku=1:nu
	bdum=[zeros(1,nk(j,ku)) b(rs+1:rs+nb(j,ku))]; rs=rs+nb(j,ku);
	yhat=yhat+filter(bdum,f,zv(:,ku));
	end
	V(1,j)=(zv(:,1)-yhat)'*(zv(:,1)-yhat)/(Ncap-nmax);
	V(2:nl+1,j)=nn(j,:)';
	if p~=0, V(nl+2,j)=log(cond(RR));,end
end
V(1,nm+1)=Ncap-nmax;


