% this file is for running regem with stepwise precedure/reconstruction of
% low-frequency band.

for istep=1:19
    
    climate=load('climate');% climate data used for first step of run, so it should have
                           % NaNs for the first decades of yrs, if no, will also be ok.
    eofnumb=load('eofnumb'); % determining the # of EOFs                      
    proxy=load(strcat('proxy',num2str(istep)));
    [ni,mi]=size(climate);
    [np,mp]=size(proxy);
    input(1:np,1:mi+mp)=NaN;
    input(1+np-ni:np,1:mi)=climate(1:ni,1:mi); % add climate data
    input(1:np,mi+1:mi+mp)=proxy(1:np,1:mp); % add proxy data
    clear proxy
    
    [n,m]=size(input);
    
    weight=[];
    X=input;

                %
                % Options for regem
                %
%                OPTIONS.regress = 'mridge';
                OPTIONS.regress = 'ttls';
                OPTIONS.regpar = eofnumb(istep)
                OPTIONS.stagtol = 5e-3;
                OPTIONS.maxit = 100;
                OPTIONS.inflation = 1;
                OPTIONS.disp = 1;
                OPTIONS.relvar_res = 0.05;
                OPTIONS.minvarfrac=0.85;

                %OPTIONS.Xcmp = keepers;
                %
%                [X, M, C, Xerr, Xmis, dXmis,peff] = newregem(X, OPTIONS, weight);
                [X, M, C, Xerr, Xmis, dXmis] = newregem(X, OPTIONS, weight);



  
    recon=X(1:np,1:mi);
    save(strcat('recon',num2str(istep)),'recon','-ASCII');
    
    nnd=[];
    save(strcat('step',num2str(istep)),'nnd','-ASCII');
    clear all
end