function [recon,rad,proxy,proxyinfo]=oie_v1_3(lat,lon,inst,preproxy,latlong)
%{
Optimal Information Extraction (OIE) version 1.3
Inputs:
  lat = latitude of the target.
  lon = longitude of the target.
  inst = the instrumental data
  preproxy = the proxy records (2000yr × n proxies)
  latlong = the infomation of all proxies includes the latitude, longitude,
  begining year, and sequence number (ID)
Outputs:
  recon = the two reconstructions for each grid: (2000yr × 3)
    col 1: the first one is the median values of the linear gression result     
    using the weight method 1 and the weight method 2;
    col 2: the secord one is the median values of the E-LOC result using the 
    weight method 1 and the weight method 2
    col 3: the number of the screened predictors per year
  rad = the search radius
  proxy = the screened predictors (2000yr × n screened proxies)
  proxyinfo = the information of the screened predictors includes the 
  latitude, longitude, distance to the target, begining year, and sequence 
  number (ID)

  Written by
  Feng SHI
  Institute of Geology and Geophysics, Chinese Acedamy of Sciences
  Beijing, China, 2017
  E-mail:shifeng@mail.iggcas.ac.cn   

  Reference:
  Shi, F., Zhao, S., Guo, Z., Goosse, H., Yin, Q.: Multi-proxy 
  reconstructions of May-September precipitation field in China over 
  the past 500 years, Climate of the Past, 2017

%}
%%
%step 1: search the proxy

yrb=1961;yre=2000;%calibration period

%calculate the correlation coeeficient between all proxies and inst
r=nan(2,2);
for i=1:size(preproxy,2)
    [r(i,1),r(i,2)]=corr(preproxy(yrb:yre,i),inst(1:yre-yrb+1,1));
    [r(i,3),r(i,4)]=corr(detrend(preproxy(yrb:yre,i)),detrend(inst(1:yre-yrb+1,1)));
end
inddd=find(r(:,2)<0.1&r(:,1)>0&r(:,4)<0.1&r(:,3)>0);
preproxy=preproxy(:,inddd);%all possible predictors
%calculate the distance between the predictors and the target    
rpoint=nan(length(inddd)+1,2);
rpoint(1,1)=lon;rpoint(1,2)=lat;
rpoint(2:end,1)=latlong(inddd,2);
rpoint(2:end,2)=latlong(inddd,1);
Dists_GC_Mat=EarthDistances(rpoint);
rpoint(:,3)=Dists_GC_Mat(:,1);
rpoint(2:end,4)=latlong(inddd,3);
rpoint(2:end,5)=1:size(rpoint,1)-1;
t1=sortrows(rpoint(2:end,:),3);
indrad= t1(:,3)>3500;%exclude the predictor with large radius
t1(indrad,:)=[];

ind_1=find(t1(:,4)<=1470);
ind_2=find(t1(:,4)==min(t1(:,4)));
%divide four conditions 
%the radius range must be no less than 3500km
%class one
%1)the number of predictor for every grid is no less than five 
%2)the time span for every grid covers AD 1470
%class two
%1)the number of predictor for every grid is no less than five 
%2)the time span for every grid can not cover AD 1470
%3)the maxmium time span is selected
%class three
%1)the number of predictor for every grid is less than five 
%2)the time span for every grid covers AD 1470
%3)the maxmium number of predictor is selected
%class four (no one)
%1)the number of predictor for every grid is less than five 
%2)the time span for every grid can not cover AD 1470
%3)the maxmium time span and the maxmium number of predictor are selected
%finally, check rad to ensure its range is [450,3500]

switch isempty(ind_1)
    case 0
    switch ind_1(1)>=5
        case 1
            indnew=t1(1:ind_1(1),5); 
            rad=t1(ind_1(1),3);%search radius 
        case 0
            switch size(t1,1)>=5
                case 1
                indnew=t1(1:5,5); 
                rad=t1(5,3);%search radius   
                case 0
                indnew=t1(1:size(t1,1),5); 
                rad=t1(end,3);%search radius 
            end
    end          
    case 1
     switch ind_2(1)>=5
        case 1
            indnew=t1(1:ind_2(1),5); 
            rad=t1(ind_2(1),3);%search radius 
        case 0
            switch size(t1,1)>=5
                case 1
                indnew=t1(1:5,5); 
                rad=t1(5,3);%search radius   
                case 0
                indnew=t1(1:size(t1,1),5); 
                rad=t1(end,3);%search radius 
            end
    end          
end

proxy=preproxy(:,indnew);%the final used proxy records
proxyinfo=t1(1:length(indnew),:);
proxyinfo(:,5)=inddd(indnew,1);
pro_num_yr=NaN(2000,1);

for i=1:2000
    proxynum1=find(proxy(i,:)>-999);
    pro_num_yr(i,1)=length(proxynum1);
end

indd=unique(proxyinfo(:,4));
proxycell = cell(length(indd),1);
for iii=1:length(indd)
    x= proxy(indd(iii,1),:)>-9999;
    proxycell{iii,1}=proxy(1:2000,x);
    clear x
end

%%
%step 2 reconstruction
n=length(indd);
recon1=NaN(yre,n);recon2=NaN(yre,n);
recon3=NaN(yre,n);recon4=NaN(yre,n);
%%%%%%input proxy inst
for istep=1:n
    yy=proxycell{istep,1};
[row1,col]=size(yy);
yyy=NaN(row1,col);
zz1=NaN(row1,col);
zz2=NaN(row1,col);
r=NaN(col,2);
%%instrumental data
y=inst(1:yre-yrb+1,1);
inststd=std(y);instmean=mean(y);

% the weighting of each record using the correlation with the instrumental data
for iiii=1:col
    [r(iiii,1),r(iiii,2)]=corr(yy(yrb:yre,iiii),y);
    yyy(:,iiii)=(yy(:,iiii)-mean(yy(yrb:yre,iiii)))/std(yy(yrb:yre,iiii))*inststd+instmean;
end
    rr(1,1)=nansum(r(:,1));
    rr(1,2)=nansum(r(:,1).^2); 
for k=1:col
    for j=1:row1      
    zz1(j,k)=yyy(j,k)*r(k,1)/rr(1,1);
    zz2(j,k)=yyy(j,k)*r(k,1).^2/rr(1,2);
    end
end
for i=1:2000
    ind=find(zz1(i,:)>-9999);
    if ind>0
z1(i,1)=nansum(zz1(i,:));
z2(i,1)=nansum(zz2(i,:));
    else
z1(i,1)=nan;
z2(i,1)=nan;     
    end
end

%clear zz1 zz2 r rr yy
x1=z1(yrb:yre,1);
x2=z2(yrb:yre,1);
p1 = polyfit(x1,y,1);p1_2 = polyfit(y,x1,1);
p2 = polyfit(x2,y,1);p2_2 = polyfit(y,x2,1);
%the linear gression result using the weight method 1
recon1(indd(istep,1):end,istep)=polyval(p1,z1(indd(istep,1):end,1));
%the linear gression result using the weight method 2
recon2(indd(istep,1):end,istep)=polyval(p2,z2(indd(istep,1):end,1));
%the E-LOC result using the weight method 1
recon3(indd(istep,1):end,istep)=0.5*(p1(1,1)+1/p1_2(1,1))*...
    z1(indd(istep,1):end,1)+0.5*(p1(1,2)-p1_2(1,2)/p1_2(1,1));
%the E-LOC result using the weight method 2
recon4(indd(istep,1):end,istep)=0.5*(p2(1,1)+1/p2_2(1,1))*...
    z2(indd(istep,1):end,1)+0.5*(p2(1,2)-p2_2(1,2)/p2_2(1,1));
end
recon=nanmedian([recon1,recon2],2);
recon(:,2)=nanmedian([recon3,recon4],2);
recon(:,3)=pro_num_yr;
 

