function [Yhat,PI]=nnvalid(method,NetDef,NN,W1,W2,par1,par2,par3)
%  NNVALID
%  ------- 
%          Validate a neural network input-output model of a dynamic system.
%          I.e., a network model which has been generated by NNARX, NNRARX,
%          NNARMAX1+2, NNRARMX1+2, or NNOE.
%
%          The following plots are produced:
%          o  Observed output together with predicted output
%          o  Prediction error
%          o  Auto-correlation function of prediction error and cross-
%             correlation between prediction error and input
%          o  A histogram showing the distribution of the prediction errors
%          o  Coefficients of extracted linear models
%
%  Call: 
%  Network generated by NNARX (or NNRARX):
%           [Yhat,NSSE] = nnvalid('nnarx',NetDef,NN,W1,W2,Y,U)
%
%  Network generated by NNARMAX1 (or NNRARMAX1):
%           [Yhat,NSSE] = nnvalid('nnarmax1',NetDef,NN,W1,W2,C,Y,U)
%
%  Network generated by NNARMAX2 (or NNRARMX2):
%           [Yhat,NSSE] = nnvalid('nnarmax2',NetDef,NN,W1,W2,Y,U)
%
%  Network generated by NNOE:
%           [Yhat,NSSE] = nnvalid('nnoe',NetDef,NN,W1,W2,Y,U)
%
%  Network generated by NNARXM:
%           [Yhat,NSSE] = nnvalid('nnarxm',NetDef,NN,W1,W2,Gamma,Y,U)
%
%  NB: For time-series, U is left out!
% 
%  Programmed by : Magnus Norgaard, IAU/IMM, Technical University of Denmark
%  LastEditDate  : June 16, 1997


% >>>>>>>>>>>>>>>>>>>>>>>>>>>>      GET PARAMETERS     <<<<<<<<<<<<<<<<<<<<<<<<<<<< 
skip = 1;
if strcmp(method,'nnarx') | strcmp(method,'nnrarx'),
  mflag=1;
  Y=par1;
  if exist('par2') U=par2; end
elseif strcmp(method,'nnarmax1') | strcmp(method,'nnrarmx1'),
  mflag=2;
  C=par1;
  Y=par2; 
  if exist('par3') U=par3; end
elseif strcmp(method,'nnarmax2') | strcmp(method,'nnrarmx2'),
  mflag=3;
  Y=par1; 
  if exist('par2') U=par2; end
elseif strcmp(method,'nnoe'),
  mflag=4;
  Y=par1;
  U=par2;
elseif strcmp(method,'nnarxm'),
  mflag=5;
  Gamma = par1;
  Y = par2;
  if exist('par3') U=par3; end
else
  disp('Unknown method!!!!!!!!');
  break
end


% >>>>>>>>>>>>>>>>>>>>>>>>>>>>     INITIALIZATIONS     <<<<<<<<<<<<<<<<<<<<<<<<<<<< 
Ndat     = length(Y);                   % # of data
na = NN(1);

% ---------- NNARX model ----------
if mflag==1 | mflag==4,
  if length(NN)==1                      % nnar model
    nb = 0;
    nk = 0;
    nu = 0;
  else                                  % nnarx or nnoe model
    [nu,N] = size(U); 
    nb = NN(2:1+nu); 
    nk = NN(2+nu:1+2*nu);
  end
  nc = 0;

% --------- NNARMAX1 model --------
elseif mflag==2 | mflag==3,
  if length(NN)==2                      % nnarma model
    nc     = NN(2);
    nb     = 0;
    nk     = 0;
    nu     = 0;
  else                                  % nnarmax model
    [nu,Ndat]= size(U); 
    nb     = NN(2:1+nu);
    nc     = NN(2+nu);
    nk     = NN(2+nu+1:2+2*nu);
  end

% ---------- NNARXM model ----------
elseif mflag==5,
  [outputs,Ndat]  = size(Y);
  [outputs,NNn]   = size(NN);
  na = NN(:,1);
  if NNn==1
    nb = 0;                          % nnar model
    nk = 0;
    nu = 0;
    nab=na;
  else
    [nu,Ndat] = size(U); 
    nb     = NN(:,2:1+nu);           % nnarx model
    nk     = NN(:,2+nu:1+2*nu);
    if nu>1,
      nab  = na + sum(nb')';
    else
      nab    = na+nb;
    end
  end
  nc = 0;
  nmax = max(max([na nb+nk-1]));
  if isempty(Gamma), Gamma=eye(outputs); end
end
  
  
% --------- Common initializations --------
if mflag>=1 & mflag<=4,
  nmax     = max([na,nb+nk-1,nc]);      % 'Oldest' signal used as input to the model
  nab      = na+sum(nb);                % na+nb
  nabc     = nab+nc;                    % na+nb+nc
  outputs     = 1;                        % Only MISO models considered
end
N        = Ndat - nmax;                 % Size of training set
L_hidden = find(NetDef(1,:)=='L')';     % Location of linear hidden neurons
H_hidden = find(NetDef(1,:)=='H')';     % Location of tanh hidden neurons
L_output = find(NetDef(2,:)=='L')';     % Location of linear output neurons
H_output = find(NetDef(2,:)=='H')';     % Location of tanh output neurons
[hidden,inputs] = size(W1);
inputs          = inputs-1;
E        = zeros(outputs,N);
y1       = zeros(hidden,N);
Yhat     = zeros(outputs,N);


% >>>>>>>>>>>>>>>>>>>>  CONSTRUCT THE REGRESSION MATRIX PHI   <<<<<<<<<<<<<<<<<<<<<
PHI = zeros(sum(nab),N);
jj  = nmax+1:Ndat;
index = 0;
for o=1:outputs,
  for k = 1:na(o), PHI(k+index,:)    = Y(o,jj-k); end
  index = index+na(o);
  for kk = 1:nu,
    for k = 1:nb(o,kk), PHI(k+index,:) = U(kk,jj-k-nk(o,kk)+1); end
    index = index + nb(o,kk);
  end
end



% >>>>>>>>>>>>>>>>>>>>>>>>>>   COMPUTE NETWORK OUTPUT   <<<<<<<<<<<<<<<<<<<<<<<<<<<
% ---------- NNARX model ----------
if mflag==1 | mflag==5,
  Y  = Y(:,nmax+1:Ndat);
  h1 = W1*[PHI;ones(1,N)];  
  y1(H_hidden,:) = pmntanh(h1(H_hidden,:));
  y1(L_hidden,:) = h1(L_hidden,:);
    
  h2 = W2*[y1;ones(1,N)];
  Yhat(H_output,:) = pmntanh(h2(H_output,:));
  Yhat(L_output,:) = h2(L_output,:);
if mflag==5, Yhat=sqrtm(Gamma)*Yhat; end
  E        = Y - Yhat;                    % Error between Y and deterministic part
  SSE      = sum(sum(E.*E));              % Sum of squared errors (SSE)
  PI       = SSE/(2*N);                   % Performance index


% --------- NNARMAX1 model --------
elseif mflag==2,
  Y  = Y(nmax+1:Ndat);
  h1 = W1*[PHI;ones(1,N)];  
  y1(H_hidden,:) = pmntanh(h1(H_hidden,:));
  y1(L_hidden,:) = h1(L_hidden,:);
    
  h2 = W2*[y1;ones(1,N)];
  Yhat(H_output,:) = pmntanh(h2(H_output,:));
  Yhat(L_output,:) = h2(L_output,:);

  Ebar     = Y - Yhat;                    % Error between Y and deterministic part
  E        = filter(1,C,Ebar);            % Prediction error
  Yhat     = Y - E;                       % One step ahead prediction

  SSE      = E*E';                        % Sum of squared errors (SSE)
  PI       = SSE/(2*N);                   % Performance index


% --------- NNARMAX2 model --------
elseif mflag==3,
  Y  = Y(nmax+1:Ndat);
  PHI_aug=[PHI;zeros(nc,N);ones(1,N)];
  y1 = [y1;ones(1,N)];
  N2=N+1-skip;
  for t=1:N,
    h1 = W1*PHI_aug(:,t);  
    y1(H_hidden,t) = pmntanh(h1(H_hidden));
    y1(L_hidden,t) = h1(L_hidden);    

    h2 = W2*y1(:,t);
    Yhat(H_output,t) = pmntanh(h2(H_output,:));
    Yhat(L_output,t) = h2(L_output,:);

    E(:,t) = Y(:,t) - Yhat(:,t);          % Prediction error
    for d=1:min(nc,N-t),
      PHI_aug(nab+d,t+d) = E(:,t);
    end
  end
  SSE      = E(skip:N)*E(skip:N)';        % Sum of squared errors (SSE)
  PI       = SSE/(2*N2);                  % Performance index


% ---------- NNOE model ----------
elseif mflag==4,
  Y  = Y(nmax+1:Ndat);
  PHI_aug=[PHI;ones(1,N)];
  y1 = [y1;ones(1,N)];
  N2=N+1-skip;
  for t=1:N,
    h1 = W1*PHI_aug(:,t);;  
    y1(H_hidden,t) = pmntanh(h1(H_hidden));
    y1(L_hidden,t) = h1(L_hidden);    

    h2 = W2*y1(:,t);
    Yhat(H_output,t) = pmntanh(h2(H_output,:));
    Yhat(L_output,t) = h2(L_output,:);

    for d=1:min(na,N-t),
      PHI_aug(d,t+d) = Yhat(:,t);
    end
  end
  E     = Y - Yhat;                       % Error between Y and deterministic part
  SSE      = E(skip:N)*E(skip:N)';        % Sum of squared errors (SSE)
  PI       = SSE/(2*N2);                  % Performance index
end


% >>>>>>>>>>>>>>>>>>>>>>>>>>      PLOT THE RESULTS      <<<<<<<<<<<<<<<<<<<<<<<<<<<
si=figure-1;


% ---------- Output, Prediction and Prediction error ----------
for ii=1:outputs,
 figure(si+ii)
 subplot(211)
 plot(Y(ii,:),'b-'); hold on
 plot(Yhat(ii,:),'r--');hold off
 xlabel('time (samples)')
 if outputs==1,
   title('Output (solid) and one-step ahead prediction (dashed)')
 else
   title(['Output (solid) and one-step ahead prediction (dashed) (output # ' ...
          num2str(ii) ')']);
 end
 grid

 subplot(212)
 plot(E(ii,:));
 title('Prediction error (y-yhat)')
 xlabel('time (samples)')
 grid
 subplot(111)
 drawnow
end

% --------- Correlation functions ----------
for ii=1:outputs,
  figure(si+outputs+ii)
  eval(['subplot(' num2str(nu+1) '11)']);
  M=min(25,N-1);
  Eauto=xcorr(E(ii,:),'coeff');
  Eauto=Eauto(N:2*N-1);
  conf=1.96/sqrt(N);
  plot([0:M],Eauto(1:M+1),'b-'); hold on
  plot([0 M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  if outputs==1
    title('Auto correlation function of prediction error')
  else
    title(['Auto correlation function of prediction error (output # ' ...
            num2str(ii) ')']);
  end
  grid
  Ecov=cov(E(ii,:));

  for i=1:nu,
    eval(['subplot(' num2str(nu+1) '1' num2str(i+1) ')']);
    Ucov=cov(U(i,1:N));
    UEcross=xcov(E(ii,:),U(i,1:N),'unbiased')/sqrt(Ecov*Ucov)';
    plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
    plot([-M M],[conf -conf;conf -conf],'r--');hold off
    xlabel('lag')
    title(['Cross correlation fct of u' num2str(i) ' and prediction error'])
    ymax=min(5*conf,max([abs(UEcross)]));
    axis([-M M -ymax ymax]);
    grid
  end
  subplot(111)
  drawnow
end

% ---------- Extract linear model from network ----------
dy2dx=zeros(outputs*(inputs+1),N);

% Matrix with partial derivative of each output with respect to each of the
% outputs from the hidden neurons
for t=1:N,
  dy2dy1 = W2(:,1:hidden);
  for j = H_output',
    dy2dy1(j,:) = W2(j,1:hidden)*(1-Yhat(j,t).*Yhat(j,t));
  end

  % Matrix with partial derivatives of the output from each hidden neurons with
  % respect to each input:
  dy1dx = W1;
  for j = H_hidden',
    dy1dx(j,:) = W1(j,:)*(1-y1(j,t).*y1(j,t));
  end

  % Matrix with partial derivative of each output with respect to each input
  dl       = (dy2dy1 * dy1dx)';
  dl(inputs+1,:)=dl(inputs+1,:)+W2(:,hidden+1)';
  dy2dx(:,t) = dl(:);
end

figure(si+2*outputs+1)
subplot(212)
plot(dy2dx(1:outputs*inputs,:)')
title('Linearized network parameters')
xlabel('time (samples)')
grid
for ii=1:outputs,
 eval(['subplot(2' num2str(outputs) num2str(ii) ')']);
 hist(E(ii,:),20)
end
eval(['subplot(2' num2str(outputs) '1)']);
title('Histogram of prediction errors')
subplot(111)
figure(si+1)