%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
% Modified PaiCO run-code for the Greenland, Arctic, Atlantic, European region
%
% Nathan Steiger, Jan 2018, Columbia University
%
%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------

disp('Initializing and loading data...')

run('../init.m')

%parpool

% LOAD PROXY DATA, BOTH THIS DATABASE AND THE SHAKUN DATABASE

opath='/home/nsteiger/MATLAB/PaiCo/lgmdeg/results/';
%mpath = fileparts(mfilename('fullpath'));
%load([mpath filesep '..' filesep 'data' filesep 'arcticatlantic.mat']);
% Remove proxies 23 (NGRIP) and 27 (GRIP) because of overlap with Shakun data
%load('arcticatlantic.mat')
%dp=data.proxy([1:22 24:26]);
%clear data

% load shakun and combine it with the Hanj dataset
%load('shakun_data_01.mat') % version includes lower and upper bound estimates
%load('shakun_data_03.mat') % THIS ONE WORKS 
%load('shakun_data_04.mat') % uses Shakun NH temps for 'instrumental' data 
%data.instrumental.noisestd=0.32; % 2 sigma unc for recon
load('shakun_data_ageEns_v1.mat'); % Updated age ensembles
data0=data; clear data

%data.proxy(end+1:(end+length(dp)))=dp; disp('Combining the two proxy datasets')
%data.proxy=data.proxy(1:2); disp('Just using Greenland proxies')
%data.proxy=data.proxy(2:end); 

% Choose the dates for the reconstruction
%data.target.times = -21000:20:2000;
%data.target.times = flip(1950:-50:-21000);
%data.target.times = flip(1950:-100:-21000);
%data0.target.times = flip(1950:-100:-25000);
data0.target.times = flip(1950:-100:-22000); % dating uncertainty blows up at about -23000
%data.target.times = flip(1950:-100:-6000);
%data.target.times = -6000:20:2000;

% Error of the instrumental time series? This value doesn't affect the mean
% but it does affect the resampled ensemble
%data.instrumental.noisestd=0.2224;
% value taken from the Arctic recons; seems to be problematic for this reconstruction



% Options for the reconstructions
options.resample.number = 0; Nb=1000; % My own bootstrap approach
%options.resample.number = 10; % original option
%options.resample.number = 10000; % add number of resamples for uncertainty quantification
%options.resample.number = 1000; % add number of resamples for uncertainty quantification

%disp(['Resampling ensembles = ' num2str(options.resample.number)])

options.resample.parallelize = true;
options.filecache = false;
options.damping = [];
options.errorTolerance = 1e-8;
options.maxIters = 1e2;
options.heuristicstart = true;
options.regcov = 100;

% COMPUTE THE RECONSTRUCTION
% (loop over this for different proxy dating ensembles?)
tracker();
tic;
% Jackknife ensemble reconstructions
reconEns=zeros(length(data0.target.times),Nb);
for i=1:Nb

   % Create a proxy structure for each jackknife ensemble
   data.instrumental=data0.instrumental;
   data.target.times=data0.target.times;
   for j=1:length(data0.proxy)
      data.proxy{j}.locations=data0.proxy{j}.locations;
      data.proxy{j}.data=data0.proxy{j}.data;
      data.proxy{j}.times=(data0.proxy{j}.ageEns(:,i))';
      data.proxy{j}.lower=(data0.proxy{j}.lowerEns(:,i))';
      data.proxy{j}.upper=(data0.proxy{j}.upperEns(:,i))';
   end

   % Randomly sample the proxies, use different age ensembles
   rsp=randsample(length(data.proxy),length(data.proxy)-1);    
   data.proxy=data.proxy(rsp);
 
   result = paico(data,options,@tracker);
   reconEns(:,i)=result.uncalibSignal;
end
toc;

% Old reconstruction option without my own bootstrap approach
%result = paico(data,options,@tracker);

% Add a time stamp to differentiate the reconstructions
timestp=strrep(datestr(now),' ','_');

% Saving data to 'results' directory
disp('Saving output...')
ofln=['arcatleurafr_r.' num2str(abs(data.target.times(1))) '_' num2str(data.target.times(end)) ...
  '_prxyN.' num2str(length(data.proxy)) '_' timestp '.mat'];
disp(['File name = ' ofln])
save([opath ofln]);









