%Set a flags here for hybrid.  
	ishybrid=1; % use 1 for hybrid 0 for nonhybrid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Load necessary data files.  The setup below is all based on .mat input
%files
    load csmprox00.mat; %Name of proxy data file.  The variable in the file is the same as the file name.
    pprox = csmprox00; %Sets the name of the proxy data file (from the previous line)
    load csmsub.mat; %Same story here for the model field
    csm = csmsub;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nyears,nproxies] = size(pprox); 
[nyears,ngrids] = size(csm);

%Set up matrix for RegEM.   
    composite=zeros(nyears,ngrids+nproxies); %First initialize a matrix of appropriate size.
	composite=composite./0;
    composite(1:1131,1:669)= csm(1:1131,1:669); %Now put in model data.
    composite(1:1131,670:669+nproxies)=pprox(1:1131,1:104);%Add the pseudo-proxies to the matrix.
	composite=composite(1:1131,1:669+nproxies); %Finally, cut composite at 1980 (for 1901-1980 calibration interval).

%Standardize everything
	means=nanmean(composite(1:1131,:));%This is for the model interval from 850(1) to 1980(1131)
	stdevs=nanstd(composite(1:1131,:));
	composite=composite-repmat(means,1131,1);%Subtract the mean.
	composite=composite./repmat(stdevs,1131,1);%Divide by the standard deviations

%For M05, standardization is over total interval.  Now erase pre-cal
%interval.
    composite(1:1006,1:669) = NaN;

%These are necessary inputs for RegEM.    
    [nyears,nvars] = size(composite);	
    weight=[];%For weight = emptyset, RegEM uses default weighting (see line 400 in regem.m)
    centerend=1006;%This is set for initial centering of the data matrix.
    %This is the end of the reconstruction period.  In this code the 
    %beginning has a default setting of 1. The data are centered from 850 
    %to 1856. See center.m for more.
    
    if ishybrid==0
        X=composite(1:nyears,:);
		% Options for regem
		OPTIONS.regress = 'mridge';
		%OPTIONS.regpar = 100;
		OPTIONS.stagtol = 5e-3;
		OPTIONS.maxit = 200;
		OPTIONS.inflation = 1;
		OPTIONS.disp = 1;
		OPTIONS.relvar_res = 0.050;
		OPTIONS.minvarfrac=0.95;
        
		%OPTIONS.Xcmp = keepers;
		
		[X, M, C, Xerr, Xmis, dXmis] = regem(X, OPTIONS, weight,centerend);
    
    elseif ishybrid==1;
        cutoff=20; %years
        %nproxies=112;
        lowlim=1./(cutoff./2)
        %start with some weights
        %tricky here, a higher number is a lower weight. The actual weight in RegEM
        %is the reciprocal of the weight here. i.e. RegEM devides by the weights.
        clear *weight*

        weight(1:nvars)=1;  %thats the pcproxies
        weight=diag(weight);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        [nrows,ncols]=size(composite)
        gotstart=zeros(1,ncols);
        gotstop=zeros(1,ncols);
        start=zeros(1,ncols);
        stop=zeros(1,ncols);
        stop(:)=nrows;
        [lowa,lowb]=butter(5,lowlim,'low');
%%%%%%%%%%%%%%
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

    for iii=1:nrows
        for jj=1:ncols
                if ~isnan(composite(iii,jj)) & gotstart(jj)~=1
                        start(jj)=iii;
                        gotstart(jj)=1;
                elseif isnan(composite(iii,jj)) & gotstart(jj)==1 & gotstop(jj)~=1
                        stop(jj)=iii-1;
                        gotstop(jj)=1;
                end
        end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Apply the zero-phase butterworth filter but only to the proxies
        highf=zeros(nrows,ncols);
        lowf=zeros(nrows,ncols);
        lowf=lowf./lowf;
        highf=highf./highf;
        
for k=1:ncols
        clear y y2
        data=composite(start(k):stop(k),k);
        lowf(start(k):stop(k),k)=lowpass(data,0.05,1,1);
        highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
end

%% Now scale lowf and highf proxies only by the percentage of variance
%each contains relative to the unfiltered series (ignore NaNs for this).
lowfpct=(nanstd(lowf).^2)./(nanstd(composite).^2);
lowfpct=1-lowfpct;
highfpct=(nanstd(highf).^2)./nanstd(composite).^2;
highfpct=1-highfpct;
highfpct(1:ncols-nproxies)=1; % keep instrumental weights=1
lowfpct(1:ncols-nproxies)=1; % keep instrumental weights 1
%The weights need to be a diagonal matrix
highweight=diag(highfpct);
highweight=highweight*weight; %combine the weights
%lowweight=diag(lowfpct);
%lowweight=lowweight*weight;  %combine the weights
%%%%%%%%%%%%%%%%%%%%%%%%
%% NOW RUN THE high FREQUENCY PART
clear lowweight lowfpct lowf k j jj jjj ii iii *pct*

X=highf(1:nyears,:);
% Options for regem
                OPTIONS.regress = 'mridge';
                %OPTIONS.regpar = 100
                OPTIONS.stagtol = 5e-2;
                OPTIONS.maxit = 200;
                OPTIONS.inflation = 1;
                OPTIONS.disp = 1;
                OPTIONS.relvar_res = 0.05;
                OPTIONS.minvarfrac=0.95;

                %OPTIONS.Xcmp = keepers;
                %
                [X, M, C, Xerr, Xmis, dXmis] = regemhigh_hpar(X,OPTIONS,highweight,centerend);
end % end the ishybrid part
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%

%Last Section
csm00_m05_high =X(1:1131,1:669); %LHS is the name of the completed Regem reconstruction.

save finished_csm00_m05_sub_high
clear C M OPTIONS X Xerr Xmis pprox echog ngrids allcenterend centerend composite constant dXmis i ishybrid mbhstepwise means nproxies nvars nyears season stdevs temp weight yearsteps 
%%%%Finished!