function V=arxstruc(z,zv,nn,maxsize);
%ARXSTRUC	Computes loss functions for families of ARX-models
%
%	V = arxstruc(Z,ZV,NN)
%
%	V: The first row of V is returned as the loss functions of the 
%	structures defined by the rows of NN. The remaining rows are 
%	returned as NN'. The last column of V contains the number of 
%	data points in Z.
%
%	Z : the output-input data Z=[y u], with y and u as column vectors.
%	Only Single-Output data are handled by this routine!
%	For a time series Z=y only.
%	ZV: the input-output data on which the validation is performed.
%	(could equal Z). The number of rows in ZV need not equal those in Z.
%
%	NN: Each row of NN is of the format [na nb nk], the orders and delays
%	of the corresponding ARX model.(See HELP ARX and HELP STRUC.)
%	For a time series, each row of NN equals na only.
%
%	See HELP SELSTRUC for analysis of V.
%	Some parameters associated with the algorithm are accessed by
%	V = arxstruc(Z,ZV,NN,maxsize)
%	See HELP AUXVAR for an explanation of maxsize.

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

[Ncap,nz]=size(z); nu=nz-1;[nm,nl]=size(nn);[Ncapv,nzv]=size(zv);
if nu>1 & nm==1,error('For just one model, use ARX instead!'),end
if nz>Ncap,error(' Data should be organized in column vectors!'),return,end
if nl~=1+2*(nz-1),disp(' Incorrect number of orders specified!'),disp('For an AR-model nn=na'),disp('For an ARX-model, nn=[na nb nk]'),
error('see above'),return,end
if nzv~=nz,error(' ZV-set should have the same number of inputs as Z!'),return,end
if nz>1, na=nn(:,1);nb=nn(:,2:1+nu);nk=nn(:,2+nu:1+2*nu);
else na=nn(:,1); nb=zeros(nm,1);nk=zeros(nm,1);end
nma=max(na);nbkm=max(nb+nk)-ones(1,nu);nkm=min(nk);
n=nma+sum((nbkm-nkm))+nu;
if nn==0,return,end
vz=1;if size(z)==size(zv),if norm(z-zv)<eps, vz=0;end,end

% *** Set up default values **
maxsdef=idmsize(Ncap,n);

if nargin<4, maxsize=maxsdef;end
if isempty(maxsize),maxsize=maxsdef;end
if maxsize<0, maxsize=maxsdef; end

% *** construct regression matrix ***
 
 nmax=max(max([na+ones(nm,1) nb+nk]'))-1;
M=floor(maxsize/n);
R=zeros(n);F=zeros(n,1);if vz,Rv=zeros(n);Fv=zeros(n,1);end
for k=nmax:M:max(Ncap,Ncapv)
	if min(Ncap,k+M)<k+1,ntz=0;else ntz=1;end
	if min(Ncapv,k+M)<k+1,ntzv=0;else ntzv=1;end
      if ntz,jj=(k+1:min(Ncap,k+M));phi=zeros(length(jj),n);end
      if vz & ntzv,jjv=(k+1:min(Ncapv,k+M));phiv=zeros(length(jjv),n);end
      for kl=1:nma, 
      if ntz,phi(:,kl)=-z(jj-kl,1);end
      if vz & ntzv ,phiv(:,kl)=-zv(jjv-kl,1);end,end
      ss=nma;
      for ku=1:nu
           for kl=nkm(ku):nbkm(ku),
	   	if ntz,phi(:,ss+kl+1-nkm(ku))=z(jj-kl,ku+1);end
	   	if vz & ntzv,phiv(:,ss+kl+1-nkm(ku))=zv(jjv-kl,ku+1);end
	   end	
           ss=ss+nbkm(ku)-nkm(ku)+1;
      end
      if ntz,R=R+phi'*phi; F=F+phi'*z(jj,1);end
      if vz & ntzv,Rv=Rv+phiv'*phiv; Fv=Fv+phiv'*zv(jjv,1);end
end
%
% *** compute loss functions ***
%
v1=zv(nmax+1:Ncapv,1)'*zv(nmax+1:Ncapv,1);
for j=1:nm
	s=[1:na(j)];rs=nma;
	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=R(s,s);
	FF=F(s);
	if vz,RRv=Rv(s,s);FFv=Fv(s);end
	TH=(RR\FF);
	if vz,
		V(1,j)=(v1-2*FFv'*TH + TH'*RRv*TH)/(Ncapv-nmax);
	else	
		V(1,j)=(v1-FF'*TH)/(Ncap-nmax);
	end
	V(2:nl+1,j)=nn(j,:)';

end;
V(1,nm+1)=Ncap-nmax;
