% RTSTRUCT Recontoolbox structure utility function
%     data = rtstruct() returns an example dataset that can be used with
%     most reconstruction methods. 
%  
%     rtstruct(data) checks if the input data conforms to the structure
%     format. Throws an error if it does not.
%     
%     Structure of data: (optional fields are marked with *)
%     data.proxy                      
%         Cell array of proxy data
%     data.proxy{i}.data              
%         Matrix of proxy record. Columns equal times and rows equal proxy 
%         records. Records with identical temporal structure may be stacked.
%         Can contain NaN's.
%     data.proxy{i}.times             
%         A row vector of mid-times for each sample. Must have equal number 
%         of columns as data. Measured in years AD.
%     *data.proxy{i}.lower            
%         Lower, or older, bound in time for each proxy sample. Measured in 
%         years AD.
%     *data.proxy{i}.upper            
%         Upper, or younger, bound in time for each proxy sample. Measured in 
%         years AD.
%     *data.proxy{i}.intargetres      
%         Boolean stating if the record is in the same resolution as the 
%         target. Some methods handle such records differently.
%     *data.proxy{i}.locations        
%         Locations of proxy records. Must have equal number of rows as data.
%         First column is degrees in latitude, second in longitude.
%     data.instrumental              
%         Struct of instrumental data.
%     data.instrumental.data 
%         Instrumental values. Columns correspond to times and rows to
%         locations. Can contain NaN's.
%     data.instrumental.times        
%         Mid-times for samples in instrumental data. Measured in years AD.
%         Currently only annual instrumental data is considered.
%     *data.instrumental.locations    
%         Locations of instrumental data. Must have equal number of rows as
%         data. First column is degrees in latitude, second in longitude.
%     data.target
%         Struct of reconstruction target.
%     data.target.times
%         Mid-times for target of reconstruction. Measured in years AD.
%         Currently only annual targets are considered.


% From the paper, section 5
%The instrumental noise standard deviation for uncertainty estimation is 
%obtained by first calculating the variance of the mean estimator for each 
%year, and then calculating the square root of the average of these variances.

% Deduce lower and upper bounds or just leave a single date?
%ulbnds='y';

% CREATE PROXY AND INSTRUMENTAL TIME SERIES MATCHING THIS FORMAT

disp('Loading proxy data...')

latlon=xlsread('shakun_2012_minimalproxydata.xlsx',2);
plat=latlon(:,1);plon=latlon(:,2);

% Just take proxies in the quandrant over the North Atlantic, Greenland, Arctic,
%  Europe, and Northern Africa
ti=find(-60 < plon & plon < 40 & plat > 20);
plat=plat(ti);plon=plon(ti);

% Load proxy data that's just in the above quadrant
D0=xlsread('EuroAtlGrnPrxyData.xlsx'); % order of the data [Temp, age, age model error 1 sigma]

% Pull out the different types of data
%dD0=D0(:,1:4:end); % depth
pD0=D0(:,2:4:end); % proxy values = Temp conversion
aD0=D0(:,3:4:end); % age years BP 
eD0=D0(:,4:4:end); % age error (1 std)

% Target depths (line between end points)
%td=linspace(dD(1),dD(end),500);

% Number of age ensembles
N=1000; disp(['Computing N = ' num2str(N) ' age ensembles...'])

data=struct();
for i=1:length(ti)

   disp(['Working on proxy ' num2str(i)])

   % Proxy locations
   data.proxy{i}.locations=[plat(i) plon(i)];

   % Pull out only legit values; safe to do because NaNs are only on the ends...
   aD=aD0(:,i);aD=aD(~isnan(aD));
   eD=eD0(:,i);eD=eD(~isnan(eD));
   pD=pD0(:,i);pD=pD(~isnan(pD));

   % Just take first 25000 years from certain proxies => dates are poor before this period
   if ismember(i,[1,2])
      mID=find(aD<=25000);
      aD=aD(mID);eD=eD(mID);pD=pD(mID);
   end

   % check that data are times are same length
   if length(pD)~=length(aD); error('NaN removal did not work...'); end

   if ismember(i,[1,2]) % Skip age ensembles for ice cores
      % Convert all ages to CE from BP
      aD=1950-aD;
   else
      % Compute the age ensembles
      [af]=makeAgeEns(aD,eD,N);
      % Convert all ages to CE from BP
      af=1950-af;
      aD=1950-aD;
   end

   % Figure out reasonable lower and upper bounds for the proxy: assume that the
   % upper and lower bounds are just the midpoints between mean dates

   % Bin high resolution data into centennial bins
   if median(abs(diff(aD)))<50
      edg0=[2000:-100:aD(end)];
      [~,aDN,bin] = histcounts(flip(aD),flip(edg0)); % determine bins; must flip for histcount
      % Remove data on end that doesn't fit in the bins
      ridx=find(bin==0);bin(ridx)=[];fpda=flip(pD);fpda(ridx)=[];
      pdaN=accumarray(bin,fpda,[],@mean)'; % mean of each bin
      pD1=flip(pdaN);
      data.proxy{i}.data=pD1; % flip time scale back
      aD1=edg0(1:end-1)-50;% mid points of edges
      data.proxy{i}.age0=aD1;    
   else % Don't bin the data
      % Save original proxy ages
      data.proxy{i}.age0=aD(:)';
      % Proxy data
      data.proxy{i}.data=pD(:)';
   end

   % Find upper and lower bounds for the age ensemble
   if median(abs(diff(aD)))<50
      lwrsEns=zeros(length(aD1),1);
      uprsEns=zeros(length(aD1),1);
      data.proxy{i}.ageEns=repmat(aD1(:),[1 N]);
   else
      lwrsEns=zeros(size(af));
      uprsEns=zeros(size(af));
      data.proxy{i}.ageEns=af;
   end
   for jj=1:N 

      % Don't worry about binned data since the time scale is so much longer
      if median(abs(diff(aD)))<50
	 dts=aD1;
      else      
	 dts=af(:,jj); 
      end
      % Pad time data for next step
      df=dts(1)+abs(dts(1)-dts(2));
      de=dts(end)-abs(dts(end)-dts(end-1));
      dp=[df dts(:)' de];

      lwrs=nan(1,length(dts));
      uprs=nan(1,length(dts));
      for ii=1:length(dts);
	 uprs(ii)=round(median(dp(ii):-1:dp(ii+1)));
	 %lwrs(ii)=median(dp(ii+1):-1:dp(ii+2))+1;
	 lwrs(ii)=round(median(dp(ii+1):-1:dp(ii+2))); % don't worry about offsetting
      end
    
      lwrsEns(:,jj)=lwrs;      
      uprsEns(:,jj)=uprs;      

   end
   
   % Save upper and lower bounds
   data.proxy{i}.lowerEns=lwrsEns;
   data.proxy{i}.upperEns=uprsEns;

end


%nlns=find(plon<0); plon(nlns)=plon(nlns)+360;

% Just take proxies in the quandrant over the North Atlantic, Greenland, Arctic,
%  Europe, and Northern Africa
%ti=find(-60 < plon & plon < 40 & plat > 20);
%data.proxy=data.proxy(ti);


%====================================================================
% Get the instrumental time series over this region
[X_obs,xlat,xlon,x_yrs] = load_obs_S('bearth',4,3);
X_obs=reshape(X_obs,length(xlat),length(xlon),length(x_yrs));

lonext=[40 (360-60)];latext=[20 85];

% Find nearest lat/lon equal to or outside bounds
ln_a=find((abs(xlon-lonext(1))-min(abs(xlon-lonext(1))))<1e-5);
ln_b=find((abs(xlon-lonext(2))-min(abs(xlon-lonext(2))))<1e-5);
% if two equal mins, pick one outside of bound
ln1=min(ln_a);ln2=max(ln_b);

% Find nearest lat/lon equal to or outside bounds
lt_a=find((abs(xlat-latext(1))-min(abs(xlat-latext(1))))<1e-5);
lt_b=find((abs(xlat-latext(2))-min(abs(xlat-latext(2))))<1e-5);
% if two equal mins, pick one outside of bound
lt1=min(lt_a);lt2=max(lt_b);

% Just take over the region
ngae=X_obs(lt1:lt2,[ln2:360 1:ln1],:);
% area weighting
A=cosd(repmat(xlat(lt1:lt2),[1 length([ln2:360 1:ln1])]));
% compute index
ngaeInd=wmean_a(ngae,A);

% estimate instrumental noise as in Hanjavari paper
n2=reshape(ngae,length(lt1:lt2)*length([ln2:360 1:ln1]),length(x_yrs));
nvar2=nanvar(n2); % variance of the field each year
data.instrumental.noisestd=sqrt(mean(nvar2)); 

% Just take instrumental data up to 2000, i = 151
data.instrumental.data=ngaeInd(1:151)';
data.instrumental.times=x_yrs(1:151);
%====================================================================


% Use Shakun NH recon for the instrumental scaling
%snh=xlsread('shakun_nh_temp.xlsx',1); % [dates data]

%data.instrumental.data=snh(:,2)';
%data.instrumental.times=snh(:,1)';

disp('Saving data...')
save('shakun_data_ageEns_v1.mat','data')

























