%function  S = pages2k_phase2_db_unpack_AMA(vers,visual)
% function pages2k_phase2_db_unpack_flat(vers)
%  INPUT:   vers (version of the database, string)
%  OUTPUT:  saved to Pages2kPhase2Database_{vers}_unpack
%
% Loads data from one of the pages2kTS.MMYYYY.mat files and extract
% data and metadata for use in climate reconstructions.
%
%  started Nov 23 2015 by Julien Emile-Geay (USC)
%  Genealogy:  same file  by J. Wang (preNov2015 directory)

%added April-March window for export (RN 2016,03,11)

vers='1_13_1';
visual='true';

% Prepare T-sensitive records
% load data and packages
addpath(genpath('../utilities'))
%load(['../../data/proxy_db/pages2kTSv' vers '.mat'])
load(['Pages2kTSv' vers '.mat'])

% Keep a record if it is interpreted as temperature AND
clim_interp = {TS.climateInterpretation_variable};
keep        = find((strcmp(clim_interp,'T')) & strcmp({TS.paleoData_useInGlobalTemperatureAnalysis}, 'TRUE'));
S     = TS(keep); nr = length(S);
not_kept    = setdiff([1:nr],keep);  % (diagnostic purposes only)

% capture interpretation direction
%temp_direc = cellfun(@lower,{S.climateInterpretation_interpDirection});

% KEEP ALL PROXY DATA!!
% disp('KEEPING ALL PROXY DATA...')
% S     = TS; nr = length(S);

% AD HOC: fix erroneous archive types
%S(78).archiveType = 'historic';
%S(81).archiveType = 'tree';
%S(270).archiveType = 'borehole';
%S(363).archiveType = 'lake sediment';

% Archive type
archive = {S.archiveType};

yearMin = nan(nr,1);
yearMax = nan(nr,1);
resMed = nan(nr,1);
resAvg = nan(nr,1);
resMax = nan(nr,1);
resMin = nan(nr,1);

% Extract time infomation first
for r = 1:nr
    t = S(r).year;
    [yearMax(r),yearMin(r)] = maxmin(floor(t));
end
year_i = min(yearMin); year_f = max(yearMax);
year   = year_i:year_f; ny = length(year);

% Aggregate proxy timseries into 3 timeseries, f(seasonal sensitivity)
proxy_ann = nan(ny,nr);
proxy_djf = nan(ny,nr);
proxy_jja = nan(ny,nr);
proxy_ama = nan(ny,nr);

% parameters for annualization
ms = 4; me = 3; % scope of annualization, from start to end month.

% parameters for SSA interpolation
op.crit_name = 'pct_var'; op.P = 1.00;
n = 0; m = 0; gappy = zeros(nr,1);gfrac=zeros(nr,1);
for r = 1:nr
    %for r = 591  %debug
    
    T = S(r);
    disp(['Importing record ' int2str(r) '/' int2str(nr)])
    
    % KLUDGE ALERT1: Skip tricky records
    if strcmp(T.dataSetName,'Aus-Duckhole.Saunders.2013')
        disp(['Skipping proxy ' num2str(r) ' because of input issues...'])
        continue 
    end

    % KLUDGE ALERT2: some records have NaNs in the chronology. Really? Yes, really.
    t     = T.year;     noNaN = ~isnan(t);
    t     = t(noNaN);   dt    = median(diff(t));
    X = T.paleoData_values(noNaN);
    
    % KLUDGE ALERT3: some records have text values stored in the
    % 'paleoData_values' field
    if iscell(X) 
        disp(['Skipping proxy ' num2str(r) ' because data values are NaNs'])
        continue
    end 
    
    % KLUDGE ALERT4: remove false annual resolution proxies
    dX0=sum(diff(X)==0); % number of values with no slope between them
    dN=sum(~isnan(X)); % number of non-nan values
    if dX0/dN > 0.1 
        % allow only 10% of values to have no slope between them. Some 
        % records start out annual, but degrade later, so we can't just
        % throw out everything that has no slope.
        disp(['Skipping proxy ' num2str(r) ' because it is false high resolution'])
        continue
    end
    
    % KLUDGE ALERT: remove proxies that are just year values!
    pcr=nancorr(t,X);
    if pcr > 0.98
        disp(['Skipping proxy ' num2str(r) ' because data was wrongly incoded'])
        continue
    end
    
    
    % KLUDGE ALERT5: some records do not encode missing values as NaN.
    % (I pity the fools!!!)
    if max(abs(X)) > 900 && median(abs(X)) < 900
        display(['Record ' num2str(r) ' may have missing values coded as -999/999'])
        X(abs(X) > 900) = nan;
    end
    %  Some records go backwards. If so, flip them so everyone is forward-looking
    if dt < 0
        t  = flipud(t(:));
        dt = -dt;
        X = flipud(X);
    end

    % KLUDGE ALERT6: check for empty records
    noNaNx = ~isnan(X);
    if sum(noNaNx)< 2
        display(['record ',num2str(r),', site ' T.dataSetName])
        %display(['record ',num2str(keep(r)),', site ' T.dataSetName])
        n = n+1;
        miss(n) = r;
    else
        % CLEAN UP AND ANNUALIZE RECORDS
        tn = t(noNaNx); Xn = X(noNaNx);
        resMed(r) = median(diff(tn));
        resAvg(r) = mean(diff(tn));
        resMax(r) = max(diff(tn));
        resMin(r) = min(diff(tn));

        % remove duplicates
        [tc,Xc,ind] = consolidator(t(:),X(:));
        % find breaks
        [~,~,segment] = find_breaks(tc,Xc);
        nseg = length(segment);
        % find first and last data points
        tmin = min(tc(segment{1})); tmax = max(tc(segment{nseg}));
        interval = find(tc >= tmin & tc <= tmax);
        ti = tc(interval); Xi = Xc(interval);
        % find breaks
        [~,~,segment] = find_breaks(ti,Xi);

        gfrac(r) = numel(Xn)/numel(X); % gap fraction

        % identify gappy records
        if nseg > 1 && resMax(r) > 10 && resMed(r) < 1 % if more than 1 segment, and large gaps
            gappy(r) = 1;
        end

        if resMed(r) < 1 % if subannual, coarse-grain to annual resolution
            if gappy(r)  % there are large gaps, proceed segment by segment
                ta = unique(floor(ti))+1;
                Xann = NaN(size(ta)); Xdjf = Xann; Xjja = Xann;Xama = Xann;
                for s = 1:nseg
                    ts = tc(segment{s}); Xs = Xc(segment{s});
                    % interpolate onto regular time grid
                    [~,~,Xr,tr] = coral_interp(ts,Xs,dt);
                    % infill using singular spectrum analysis
                    [spec,eig_vec,PC,RC,Xi,modes] = hepta_ssam(Xr,op);
                    % intra annual averages
                    [ann,ty,~]  = intra_annual_avg(Xi,tr,ms,me);
                    Xann(ismember(ta,ty)) = ann;
                    djf         = intra_annual_avg(Xi,tr,12,2);
                    Xdjf(ismember(ta,ty)) = djf;
                    jja         = intra_annual_avg(Xi,tr, 6,8);
                    Xjja(ismember(ta,ty)) = jja;
                    ama         = intra_annual_avg(Xi,tr, 4,3);
                    Xama(ismember(ta,ty)) = ama;
                end
            else
                % interpolate onto regular time grid
                mRes = round(resMed(r)*12);
                [~,~,Xr,tr] = coral_interp(tc,Xc,1,mRes);
                % infill using singular spectrum analysis
                if sum(isnan(Xr))>0
                    [spec,eig_vec,PC,RC,Xi,modes] = hepta_ssam(Xr,op);
                else
                    Xi = Xr;
                end
                % intra annual averages
                [Xann,ta,~]  = intra_annual_avg(Xi,tr,ms,me);
                Xdjf         = intra_annual_avg(Xi,tr,12,2);
                Xjja         = intra_annual_avg(Xi,tr, 6,8);
                Xama         = intra_annual_avg(Xi,tr, 4,3);

                if sum(~ismember(ta,unique(floor(tc))))~=0
                    % make sure corals are not annualized where there were original NaNs
                    Xann(~ismember(ta,unique(floor(tc)))) = nan;
                    Xdjf(~ismember(ta,unique(floor(tc)))) = nan;
                    Xjja(~ismember(ta,unique(floor(tc)))) = nan;
                    Xama(~ismember(ta,unique(floor(tc)))) = nan;
                end
            end

            %sanity check
            if visual
                fig('Annualization'), clf
                plot(tc,Xc,'marker','o','color',rgb('Gray'),'MarkerSize',8,'linestyle','none'); hold on;
                %plot(ta,Xann,'k-','linewidth',2); % ORIGINAL
                plot(ta,Xama,'k-','linewidth',2);
                if resMed(r) < 1
                    plot(ta,Xdjf,'b-',ta,Xjja,'r-');
                end
                hold off; lab{1} = 'Raw'; lab{2} = 'Annualized'; legend(lab{:}), legend boxoff
                ylab = [T.paleoData_variableName ' (' removeLeadingAndTrailingSpaces(T.paleoData_units) ')'];
                site_n = strrep(T.dataSetName,'_','\_');
                ttl = [int2str(r),') ', archive{r},', ', site_n];
                fancyplot_deco(ttl,'Year (CE)',ylab,14);
                filen=['../figs/annualize/pages_2k_phase2_record_', sprintf('%03d',r) '.pdf'];
                %export_fig(filen,'-r100','-cmyk','-painters','-nocrop') % NJS: ISN'T WORKING ON MY MACHINE...
                print(filen,'-dpdf','-painters')
            end
        elseif resMed(r)>=1 && resMin(r) < 1
            ta = [min(floor(tc)):max(floor(tc))]';
            F = griddedInterpolant(tc,Xc,'linear'); Xann = F(ta);
            Xdjf = Xann; Xjja = Xann;Xama = Xann;
        else
            [tc,Xc,ind] = consolidator(tn(:),Xn(:));
            ta = round(tc);
            Xann = Xc; Xdjf = Xc; Xjja = Xc; Xama = Xc;
        end
        if strcmp(archive{r},'ice core')
            proxy_ann(ismember(year,ta),r) = Xann;  % As per Eric Steig: diffusion renders subannual averages meaningless
            proxy_djf(ismember(year,ta),r) = Xann;
            proxy_jja(ismember(year,ta),r) = Xann;
            proxy_ama(ismember(year,ta),r) = Xama;
        else
            proxy_ann(ismember(year,ta),r) = Xann;
            proxy_djf(ismember(year,ta),r) = Xdjf;
            proxy_jja(ismember(year,ta),r) = Xjja;
            proxy_ama(ismember(year,ta),r) = Xama;
        end
        clear Xann Xdjf Xjja Xama
    end
end

% special case of PacificGBRWei2009: resolution too darn variable to handle
% Wei, G., M.T. McCulloch, G. Mortimer, W. Deng, and L. Xie. 2009. Evidence for ocean acidification in the Great Barrier Reef of Australia. Geochimica et Cosmochimica Acta, vol. 73, pp. 2332-2346. doi:10.1016/j.gca.2009.02.009
names = {S.dataSetName};
iWei = find(strcmpi(names,'Ocean2kHR-PacificGBRWei2009')); %2 hits
for r = iWei
    T = S(r);
    tc = T.year;
    Xc = T.paleoData_values;
    resMed(r) = 5;
    resAvg(r) = 5;
    resMax(r) = 5;
    resMin(r) = 5;
    ta = [min(tc):5:max(tc)]; % 5y
    Xm = nmean(Xc);
    X5 = pchip(tc,Xc-Xm,ta)+Xm;
    ty = [min(tc):1:max(tc)]; % yearly
    Xann = pchip(ta,X5,ty);

    % save to proxy matrices
    proxy_ann(ismember(year,ty),r) = Xann;
    proxy_djf(ismember(year,ty),r) = Xann;
    proxy_jja(ismember(year,ty),r) = Xann;
    proxy_ama(ismember(year,ty),r) = Xann;
end
% restrict to Common Era
yearCE = year>0;
year  = year(yearCE);
proxy_ann = proxy_ann(yearCE,:);
proxy_djf = proxy_djf(yearCE,:);
proxy_jja = proxy_jja(yearCE,:);
proxy_ama = proxy_ama(yearCE,:);

% Coordinates
p_lat = [S.geo_meanLat];
p_lon = [S.geo_meanLon];

%-----------------------------
% FURTHER PROXY DATA CLEAN UP
%-----------------------------

% Remove proxies that are just NaNs because they were skipped above
ni = find(sum(isnan(proxy_ama),1)==size(proxy_ama,1));
archive(ni)=[];p_lat(ni)=[];p_lon(ni)=[];proxy_ann(:,ni)=[];
proxy_ama(:,ni)=[];resMin(ni)=[];

% Remove longer than annual proxies
li=find(resMin>1);
archive(li)=[];p_lat(li)=[];p_lon(li)=[];proxy_ann(:,li)=[];
proxy_ama(:,li)=[];

% Remove proxy repeats (very highly correlated proxies)
rtol=0.98;
[proxy_ama,p_lat,p_lon,archive] = rmv_repeats(proxy_ama,p_lat,p_lon,archive,rtol);


clear TS D tr Xr tc Xc ts Xs spec eig_vec PC RC  Xi modes Xn tn n m r Xann

% file name corresponding to current time stamp
timestp=datestr(now);
timestp=strrep(timestp,' ','_'); % remove space

%save(['pages2kTSv' vers '_unpack_AMA_allproxies_njs_' timestp '.mat'],...
%    'archive','p_lat','p_lon','proxy_ama','year')
save(['pages2kTSv' vers '_unpack_AMA_Trecon_njs_' timestp '.mat'],...
    'archive','p_lat','p_lon','proxy_ama','year')
%save(['../../data/proxy_db/pages2kTSv' vers '_unpack_AMA.mat'])
%save(['pages2kTSv' vers '_unpack_AMA.mat'])





