HIGHR2.DOG

function [ikey,isurr,rsurr,nsurr,errflg]=highr2(R,j1,j2,j3,mlap,rmin)
% highr2: highest correlation between tree-ring series; surrogate series
% CALL: [ikey,isurr,rsurr,nsurr,errflg]=highr2(R,j1,j2,j3,mlap,rmin);
%
% Meko 4-19-97
%
% Subfunction written for use with surrgt?.m functions that
% identify surrogate tree-ring series
%
%******************* IN *********
%
% R (mR x 4)r correlation, sample size, series 1 index, series 2 index
%    as obtained from r4sov.m
% j1 (1 x nsers)L  pointer to completed series
% j2 (1 x nsers)L  pointer to yet to be completed series
% j3 (1 x nsers)L  pointer to special (specified series)
% mlap (1 x 1)i  minimum number of acceptable overlap years for
%    a correlation coefficient
% rmin (1 x 1)r minimum acceptable correlation
%
%j
%****************** OUT *********************
%
% ikey (1 x 1)i  index of newly completed series
% isurr (1 x 1)i index of its surrogate
% rsurr (1 x 1)r  correlation coefficient between key and surrogate
% nsurr (1 x 1)i  sample size (# obs) for rsurr
% errflg (1 x 1)i error flag
%    0 = no problem
%    1 = no pairs of j1,j2 series in cols 3,4 of R, after
%      ruling out j3 series
%    2 suitable pairs exist, but none with at least mlap observations
%      overlap for the computed r
%    3 suitable pairs exist, and overlap OK for at least 1, but
%      none have high enough correl coef, as specified by rmin
%
%
%***************** NOTES **********************************
%
% The correlations and associated info, R, are assumed to 
% have been generated by r4sov.m.  R is therefore assumed not
% to have any NaN elements.  highr2.m checks anyway and reports an
% error if R contains any NaNs
%
% Typically, the j1,j2,j3 logical series correspond to pointers
% to time series:
% j1 -- 'completed' series.  Those either covering whole 
%   calibration period, or if not, those previously filled in
%   for the calibration period 
% j2 -- 'waiting' series.  Those not yet with surrogate data
%    generated for the calibration period
% j3 -- 'specified' series.  Series that the user does not want
%   correlations to determine surrogate for.  The user instead
%   has specified the surrogates based on other knowledge.

errflg=[]; %initialize error flag
rsurr=[];
nsurr=[];
isurr=[];
ikey=[];

% Check inputs
[mR,nR]=size(R);
if nR~=4,
   error('R must be 4-col matrix');
end
if any(any(isnan(R)));
   disp('R, as generated by r4sov.m, should not contain any NaNs')
   error('NaN elements in R');
end



Ltemp=[size(j1,1)==1 size(j2,1)==1 size(j3,1)==1];
if ~all(Ltemp);
   error('j1, j2, j3 must all be row vectors');
end
ntemp=size(j1,2);
Ltemp=[ntemp==size(j2,2) ntemp==size(j3,2)];
if ~all(Ltemp);
   error('j1,j2,j3 not all same col size');
end
Ltemp=[islogical(j1) islogical(j2) islogical(j3)];
if ~all(Ltemp);
   error('j1,j2,j3 not all logical');
end

% Check that j1, j2 mutually exclusive
Ltemp=j1 & j2;
if any(Ltemp),
   error('Some element is turned on in both j1 and j2');
end

% Check that have at least one j2 series and at least one
% j1 series, and that those series are not j3 series
Ltemp=j1 & ~j3;
if ~any(Ltemp);
   error('No j1 series, or none that are not also j3');
end
Ltemp=j2 & ~j3;
if ~any(Ltemp);
   error('No j2 series or none that are not also j3');
end



%**************************************************************
% Identify rows of R corresponding to pairs of series such
% that (1) one series is a j1 series but not a j3 series, and
% (2) the other series is a j2 series but not a j3 series, (3)
% the sample size for correlation is at least mlap observations,
% and (4) the correlation is at least rmin

% Acceptable sample size?
L1=R(:,2)>=mlap;

% Accepable minimum r
L2=R(:,1)>=rmin;

% Store cols 3 and 4 of R in a and b
a=R(:,3); % cv of indices, first of pair
b=R(:,4); % cv of indices, second of pair


% First or second series (col 3 or 4 of R) not a j3 series
if ~any(j3);
   L3a=logical(zeros(mR,1));
   L3b=L3a;
else
   numj3=sum(j3);
   j3f=find(j3); % a rv of indices of j3 series
   j3f=repmat(j3f,mR,1); % expand to a matrix, same row-size as R
   
   % Handle col 3 of R
   A=repmat(a,1,numj3); % expand to same size as j3f
   Ltemp=A==j3f;
   if numj3>1;
      L3a=(any(Ltemp'))';
   else
      L3a=Ltemp; % special case for vector instead of matrix
   end
   
   
   % Handle col 4 of R
   B=repmat(b,1,numj3); % expand to same size as j3f
   Ltemp=B==j3f;
   if numj3>1
      L3b=(any(Ltemp'))';
   else
      L3b=Ltemp;
   end
end


% Col 3 of R is a j1 series
numj1=sum(j1);
j1f=find(j1);
j1f=repmat(j1f,mR,1);

A=repmat(a,1,numj1);
Ltemp=A==j1f;
if numj1>1;
   L4a=(any(Ltemp'))';
else
   L4a=Ltemp;
end


% Col 4 of R is a j1 series
B=repmat(b,1,numj1);
Ltemp=B==j1f;
if numj1>1;
   L4b=(any(Ltemp'))';
else
   L4b=Ltemp;
end


% Col 3 of R is a j2 series
numj2=sum(j2);
j2f=find(j2);
j2f=repmat(j2f,mR,1);

A=repmat(a,1,numj2);
Ltemp=A==j2f;
if numj2>1;
   L5a=(any(Ltemp'))';
else
   L5a=Ltemp;
end


% Col 4 of R is a j2 series
B=repmat(b,1,numj2);
Ltemp=B==j2f;
if numj2>1;
   L5b=(any(Ltemp'))';
else
   L5b=Ltemp;
end


% Now we have these pointers to rows of R:
% L1 -- at least mlap observations for computation of R
% L2 -- r at least as large as rmin
% L3a -- j3 series in col 3
% L3b -- j3 series in col 4
% L4a  -- j1 series in col 3
% L4b -- j1 series in col 4
% L5a -- j2 series in col 3
% L5b -- j2 series in col 4


% One a j1 series, one a j2 series, and neither a j3 series
L6=((L4a | L4b) & (L5a | L5b)) & ~L3a & ~L3b;
if sum(L6)==0;
   disp('No j1,j2 pairs in R after excluding j3 series');
   pause(2);
   errflg=1;
   return
end

% Constraint L6 and: minimum overlap of mlap observations
L7=L6 & L1;
if sum(L7)==0;
   disp('Acceptable j1,j2 pairs, but not enought overlap in yr');
   pause(2);
   errflg=2;
   return
end

% Constraint L7 and: correlation at least rmin
L8=L7 & L2;
if sum(L2)==0;
   disp('Bombed out because correlation coef lower than rmin');
   pause(2)
   errflg=3
   return
end


%***************************************************
%
% L8 is the pointer to rows of R for which to find the max correl

f8=find(L8);  % row index in R of candidate rows

r1=R(L8,1); % cv of correlation coefficients
[rmax, imax]=max(r1);  % highest correlation coef, and row index to r1

irow = f8(imax); % row of R containing highest correlation
rsurr=R(irow,1);
nsurr=R(irow,2);
var3=R(irow,3);
var4=R(irow,4);

% Know that var3 and var4 are a j1,j2 pair, but do not know
% which is the j1 and which is the j2
if any(find(j1)==var3);
   ikey=var4;
   isurr=var3;
else
   ikey=var3;
   isurr=var4;
end
errflg=0;