function [yh,fit] = compare(z,th,m,nr,za)
% COMPARE  Compares the simulated/predicted output with the measured output
%
%	YH = compare(Z,TH,M)
%
%	Z : The output - input data for which the comparison is made (the
%	    validation data set.
%	TH: The model in the THETA format (see help theta)
%	M : The prediction horizon. Old outputs up to time t-M are used to
%	    predict the output at time t. All relevant inputs are used.
%	    M = inf gives a pure simulation of the system.(Default M=inf).
%	YH: The resulting simulated/predicted output.
%	COMPARE also plots YH together with the measured output in Z, and
%	    displays the mean square fit between these two signals.
%	    (solid/red is YH. dashed/green is measured output)
%
%	[YH,FIT] = compare(Z,TH,M,SAMPNR,LEVELADJUST)  
%	gives access to some options:
%	FIT: The mean square fit.
%	SAMPNR: The sample numbers from Z to be plotted and used for the
%	   computation of FIT. (Default: SAMPNR = all rows of Z)
%	LEVELADJUST: 'yes' adjusts the first values of YH and Z(:,1) to
%	  zero before plot and computation of FIT. 'no' is default.

%	L. Ljung 10-1-89,11-02-91
%	Copyright (c) 1989-92 by the MathWorks, Inc.
%	All Rights Reserved.
	
[nrow,nc]=size(z);nu=th(1,3);ny=nc-nu;
if nargin<5,za=[];end
if nargin<4,nr=[];end
if nargin<3, m=[];end
if isempty(za),za='n';end
if isempty(nr),nr=1:nrow;end
if isempty(m),m=inf;end
nu=th(1,3); if nu==0 & m==inf
error(['For a time series, there must be a finite prediction horizon' ...
'\n(The argument M must be less than inf)']),end
if m==inf,yh=idsim(z(:,ny+1:nc),th);else yh=predict(z,th,m);end
[Ncap,ny]=size(yh);
if za(1)=='y' | za(1)=='Y'
yh(nr,:)=yh(nr,:)-ones(length(nr),1)*yh(nr(1),:);
z(nr,1:ny)=z(nr,1:ny)-ones(length(nr),1)*z(nr(1),1:ny);end
fit=[];
for ky=1:ny
plot([yh(nr,ky) z(nr,ky)])
fittemp = norm(yh(nr,ky)-z(nr,ky))/sqrt(length(nr));	
title(['Output # ',int2str(ky),' Fit: ', num2str(fittemp)])
xlabel('Red/solid: Model output,  Green/dashed: Measured output')
fit=[fit,fittemp];
if ky<ny,pause,end
end


