function [ptd_pred,Up,Vp,Yp]=do_cca_prediction(REG,ptr,ptd,ptr_anom_field,isX)

% Function to preform prediction based on precomputed CCA
% USAGE: [ptd_pred,Up,Vp,Yp]=do_cca_prediction(REG,ptr,ptd,ptr_anom_field,isX)

% Inputs: (see example below for details)
%   REG: structure containing information about CCA
%     REG.A,REG.B: are matrices of canonical coefficients
%     REG.regstats: structure of regression statistics of the canonical variables
%   ptr: structure containing information about predictor variables
%     This routine uses ptr.X (the PCs of the multivariate model state)
%   ptd: structure containing information about predictand variables
%   ptr_anom_field: array containing the predictor raw anomaly fields OR
%                   the multivariate PCs (X) of the predictor variables.
%   isX: flag determining whether ptr_anom_field is the anomaly field 
%        (isX=0 or missing) OR X (isX=1)
% Outputs:
%   ptd_pred: predicted fields

% Samar Khatiwala (spk@ldeo.columbia.edu)

ptr_numvars=length(ptr.vars);
ptd_numvars=length(ptd.vars);

if nargin<5
  isX=0;
end

if isX % ptr_anom_field is X
  X=ptr_anom_field;
  nt=size(X,1);
else % ptr_anom_field is the raw anomaly field(s)
  disp('Assembling multivariate model state') 
  nt=size(ptr_anom_field,3);
% Construct model state composed of standarized PCs of different 
% variables
  MZpc=[];
  for iv=1:ptr_numvars
    neofs=ptr.neofsv(iv);
%     PCv=eof_project(ptr_anom_field{iv},ptr.eofs{iv});
%     ZPCv=PCv*diag(1./ptr.stdev(1:neofs));
    if isfield(ptr,'pcafiles')
      if ~isempty(ptr.pcafiles{iv})
        varFile=ptr.pcafiles{iv};
      else
        varFile=[ptr.vars{iv} '_pca'];
      end
    else
      varFile=[ptr.vars{iv} '_pca'];
    end
    S=load(varFile,'EOF','SD');
    inn=find(isnan(S.EOF(:,:,1)));
    for it=1:nt % take care of land
      anom=ptr_anom_field{iv}(:,:,it);      
      anom(inn)=NaN;
      ptr_anom_field{iv}(:,:,it)=anom;
    end
    PCv=eof_project(ptr_anom_field{iv},S.EOF(:,:,1:neofs),2);
    ZPCv=PCv*diag(1./S.SD(1:neofs));
    MZpc=[MZpc ZPCv]; % (t,e)
  end
% project onto MEOFs
  if ptr_numvars==1
    X=MZpc(:,1:ptr.meof_neofs);
  else    
% project onto (truncated set of) multivariate EOFs
  X=eof_project(MZpc',ptr.meofs_eofs(:,:,1:ptr.meof_neofs),1); % (nt,ptr.meof_neofs)
end

for iv=1:ptd_numvars
  nc=size(REG.A{iv},2);
% Project onto CCA modes
  Up{iv}=X*REG.A{iv}; % (nt,nc=ptd.neofs{iv})
% predict V
  Vp{iv}=repmat(0,[nt nc]);
  for ic=1:nc
      a=REG.regstats{ic,iv}.beta(2);
      b=REG.regstats{ic,iv}.beta(1);
      Vp{iv}(:,ic)=a*Up{iv}(:,ic)+b;
  end  
% project back to qflux PCA space
  Yp{iv}=Vp{iv}*inv(REG.B{iv});
% predict qflux
% Note: it is assumed nc=neofsq
  PC=Yp{iv}*diag(ptd.stdev{iv});
  ptd_pred{iv}=zeros([size(ptd.eofs{iv},1) size(ptd.eofs{iv},2) nt]);
  for it=1:nt
    for ie=1:ptd.neofsv(iv)
      ptd_pred{iv}(:,:,it)=ptd_pred{iv}(:,:,it)+PC(it,ie)*ptd.eofs{iv}(:,:,ie);
    end
  end
end