function Yhat=kpredict(method,NetDef,NN,p,W1,W2,par1,par2,par3)
%  KPREDICT
%  --------
%          k-step ahead predictions determined by simulation of the
%          one-step ahead neural network predictor. For NNARMAX
%          models the residuals are set to zero when calculating the
%          predictions. The predictions are compared to the observed output.
%
%  Call: 
%  Network trained with NNARX, (or NNRARX):
%           Yk = kpredict('nnarx',NetDef,NN,k,W1,W2,Y,U)
%  (Likewise for NNOE and NNARMAX1+2)
%
%  Inputs:
%        See NNVALID
%        k:   Prediction horizon
%
%  Output:
%        Ypred: k-step ahead predictions.
%
%  NB! Does not work for models generated by NNIOL or NNSSIF.
% 
%  Programmed by : Magnus Norgaard, IAU/IMM, Technical University of Denmark
%  LastEditDate  : Oct. 17, 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;
else
  disp('Unknown method!!!!!!!!');
  break
end



% >>>>>>>>>>>>>>>>>>>>>>>>>>>>     INITIALIZATIONS     <<<<<<<<<<<<<<<<<<<<<<<<<<<< 
[ny,Ndat] = size(Y);                     % # of outputs and # of data
[nu,Ndat] = size(U);                     % # of inputs 
na        = NN(1);


% ---------- NNARX/NNOE model ----------
if mflag==1 | mflag==4,
  nb = NN(2:1+nu);
  nc = 0;
  nk = NN(2+nu:1+2*nu);

% --------- NNARMAX1 model --------
elseif mflag==2, 
  nb = NN(2:1+nu);
  nc     = 0;
  nk     = NN(2+nu+1:2+2*nu);

% --------- NNARMAX2 model --------
elseif mflag==3, 
  nb     = NN(2:1+nu);
  nc     = NN(2+nu);
  nk     = NN(2+nu+1:2+2*nu);
end
nmax     = max([na,nb+nk-1]);           % 'Oldest' signal used as input to the model
N        = Ndat - nmax;                 % Size of training set
nab      = na+sum(nb);                  % na+nb
nabc     = nab+nc;                      % na+nb+nc
outputs     = 1;                        % Only MISO models considered


% --------- Common initializations --------
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;
y1       = [zeros(hidden,1);1];
y2       = zeros(outputs,1);
Yhat     = zeros(outputs,N);


% >>>>>>>>>>>>>>>>>>>>>>>>>>   COMPUTE NETWORK OUTPUT   <<<<<<<<<<<<<<<<<<<<<<<<<<<
% ---------- NNARX/NNARMAX1/NNARMAX2 model ----------
if mflag==1 | mflag==2 | mflag==3,

  % -----  CONSTRUCT THE REGRESSION MATRIX PHI   -----
  PHI_aug = [zeros(nab,N);ones(1,N)];
  jj  = nmax+1:Ndat;
  for k = 1:na, PHI_aug(k,:)    = Y(jj-k); end
  index = na;
  for kk = 1:nu,
    for k = 1:nb(kk), PHI_aug(k+index,:) = U(kk,jj-k-nk(kk)+1); end
    index = index + nb(kk);
  end
  
  % -----  CONSTRUCT PHI FOR NNARMAX2  -----
  if mflag==3,
    PHI_aug=[PHI_aug(1:nab,:);zeros(nc,N);ones(1,N)];
    N2=N+1-skip;
    for t=1:N,
      h1 = W1*PHI_aug(:,t);  
      y1(H_hidden) = pmntanh(h1(H_hidden));
      y1(L_hidden) = h1(L_hidden);    

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

      E = Y(:,nmax+t) - Yhat(:,t);          % Prediction error
      for d=1:min(nc,N-t),
        PHI_aug(nab+d,t+d) = E;
      end
    end
  
  % -----  CONSTRUCT PHI FOR NNARMAX1  -----
  elseif mflag==2,
    lc=length(C)-1;
    yy1 = zeros(hidden,N); 
    h1 = W1*PHI_aug;
    yy1(H_hidden,:) = pmntanh(h1(H_hidden,:));
    yy1(L_hidden,:) = h1(L_hidden,:);

    h2 = W2*[yy1;ones(1,N)];
    Yhat(H_output,:) = pmntanh(h2(H_output,:));
    Yhat(L_output,:) = h2(L_output,:);

    Ebar     = Y(nmax+1:Ndat) - Yhat;     % Error between Y and deterministic part
    E        = filter(1,C,Ebar);          % Prediction error
    Yhat     = Y(nmax+1:Ndat) - E;        % One step ahead prediction
    PHI2=zeros(lc,N);
    for d=1:lc,
      PHI2(d,d+1:N)=E(1:N-d);
    end
  end  
  
  % ----- DETERMINE K-STEP AHEAD PREDICTION -----
  for t=1:N-p+1,
    phi=PHI_aug(1:nabc+1,t);
    if mflag==2,
      phi2=PHI2(:,t);
    end
    for k=0:p-1,
      phi(na+1:nab)=PHI_aug(na+1:nab,t+k);
      h1 = W1*phi;  
      y1(H_hidden) = pmntanh(h1(H_hidden));
      y1(L_hidden) = h1(L_hidden);    

      h2 = W2*y1;
      y2(H_output) = pmntanh(h2(H_output,:));
      y2(L_output) = h2(L_output,:);
      if mflag==2,
        y2   = y2+sum(C(2:lc+1)*phi2);
        phi2(2:lc) = phi(1:lc-1);
        phi2(1)    = 0;
      elseif mflag==3,
        phi(nab+2:nabc)=phi(nab+1:nabc-1);
        phi(nab+1)=0;
      end
      phi(2:na)=phi(1:na-1);
      phi(1)=y2;
    end
    Yhat(:,t+p-1) = y2;
  end



% ---------- NNOE model ----------
elseif mflag==4,
   Yhat = nnsimul('nnoe',NetDef,NN,W1,W2,Y,U);
   close(gcf);
end


% >>>>>>>>>>>>>>>>>>>>>>>>>>      PLOT THE RESULTS      <<<<<<<<<<<<<<<<<<<<<<<<<<<
si = figure-1;
Y = Y(:,nmax+p:Ndat);
Yhat=Yhat(:,p:N);
for k=1:outputs,
  if outputs>1,
    figure(si+k);
  end
  plot(p:N,Y(k,:),'b-'); hold on
  plot(p:N,Yhat(k,:),'r--');hold off
  xlabel('time (samples)')
  if outputs==1,
    title(['Output (solid) and ' num2str(p) '-step ahead prediction (dashed)'])
  else
    title(['Output (solid) and ' num2str(p) '-step ahead prediction (dashed)  #' int2str(k)])
  end
  grid
end