function [T1,T2,M1,M2,eVar,eVarDist,MSE,vsl_trw,vsl_cr,eVarUV] = tree_vslite_calib(...
    calPeriod,prx_i,p_lat,p_lon,prxy_yrs,pNum)
%TREE_VSLITE_CALIB Calibrate the VS-lite model with observational data
%   Compute the parameters of the VS-lite model given a calibration period
%   tree ring data, and tree ring locations. The output is vectors of
%   parameters for all the tree locations.
%
%   Nathan Steiger, Sept 2016, LDEO

addpath('/home/nsteiger/MATLAB/regem-imputation')
addpath('/home/nsteiger/MATLAB/annual-da')

%====================================
% LOAD HISTORICAL OBSERVATIONAL DATA
%====================================

disp('Loading observational data for VS-lite calibration...')

gcm_syr=1901;

mon_avg_o=[];% leave data at monthly resolution
mon_avg_f=[];
obs_dtp1='cru_ts3_tmp';% temperature
[X_t_obs,xlat,xlon] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);

obs_dtp1='cru_ts3_pre';% precipitation
[X_p_obs] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);


% Get the nearest locations to the proxies
mask=ones(size(X_t_obs(:,end))); % mask out ocean lat/lon values
mask(isnan(X_t_obs(:,end)))=NaN;
[l_p] = nearest_latlon(xlat,xlon,p_lat,p_lon,mask);

% % CALIBRATE ONLY OVER 1911 TO 1995
c_yrs=calPeriod(1):calPeriod(2);
% [~,ig] = intersect(x_yrs,c_yrs);
[~,ip] = intersect(prxy_yrs,c_yrs);


% Pull out tree ring data over calibration years
TRW=prx_i(ip,:);


% % Reshape temp and precip to feed into the tree ring model
% tv_mon=zeros(length(xlat)*length(xlon),12*length(x_yrs));
% pv_mon=zeros(length(xlat)*length(xlon),12*length(x_yrs));
% for n=1:length(x_yrs)
%     tv_mon(:,n)=reshape(X_t_obs(:,n),length(xlat)*length(xlon),1);
%     pv_mon(:,n)=reshape(X_p_obs(:,n),length(xlat)*length(xlon),1);
% end

% Just take data at proxy locations
tp_mon=X_t_obs(l_p,:);
pp_mon=X_p_obs(l_p,:);

% infill missing values with regem
if sum(isnan(tp_mon(:))>0)
    options.maxit=5; options.regress='ttls';
    tp_mon=regem(tp_mon,options);
end
if sum(isnan(pp_mon(:))>0)
    options.maxit=5; options.regress='ttls';
    pp_mon=regem(pp_mon,options);
    pp_mon(pp_mon<0)=0.01; % make sure there aren't any zeros
end

% Take only model data corresponding to calibration period
sm=12*(calPeriod(1)-gcm_syr)+1;
em=12*(calPeriod(2)-gcm_syr+1);
calYrs=length(sm:em)/12;

% Take it at this point b/c full tp_mon/pp_mon used in model
tp_mon_cal=tp_mon(:,sm:em);
pp_mon_cal=pp_mon(:,sm:em);

% Reshape calibration data for calibrating TRW ensemble
k=1;l=12;
t_cal=zeros(pNum,12,calYrs);
p_cal=zeros(pNum,12,calYrs);
for n=1:calYrs
    t_cal(:,:,n)=tp_mon_cal(:,k:l);
    p_cal(:,:,n)=pp_mon_cal(:,k:l);
    k=k+12; l=l+12;
end

% % OLD METHOD OF BIAS CORRECTING THE MODEL/OBS DATA
% % Bias-correct model precip and temperature data
% for jj=1:length(t_lon)
%     if t_lon(jj) >= 180
%         t_lon(jj)=t_lon(jj)-360;
%     end
% end
% t_seas=zeros(pNum,12);
% p_seas=zeros(pNum,12);
% parfor p=1:pNum
%     [t_seas(p,:),p_seas(p,:)]=t_pcp_wc(t_lat(p),t_lon(p));
% end
%
% % Remove mean and add in instrumental-based mean
% p_cal_d=p_cal-repmat(mean(p_cal,3),[1,1,calYrs]);
% p_bcal=p_cal_d+repmat(p_seas,[1,1,calYrs]);
%
% t_cal_d=t_cal-repmat(mean(t_cal,3),[1,1,calYrs]);
% t_bcal=t_cal_d+repmat(t_seas,[1,1,calYrs]);
%
% % correct for pcp over-correction, min pcp = 0
% p_bcal(p_bcal<0)=0;


% Use only good non-nan years in the parameter estimations
gyrs=cell(pNum,1);
for i=1:pNum; gyrs{i}=find(~isnan(TRW(:,i))); end

disp('Estimating VS-lite parameters...')

% Find calibration parameters for all proxy locations
% growth parameters, modifying lower-bound temperature growth to 7.5 degC
fp=[9 5 0 7.5; 3.5 3.5 10 24; 1.5 2.8 0.0 0.1; 1.5 2.5 0.1 0.5];
T1=zeros(pNum,1);T2=zeros(pNum,1);
M1=zeros(pNum,1);M2=zeros(pNum,1);
eVar=zeros(pNum,1); % model error variance
eVarDist=cell(pNum,1);
parfor q=1:pNum
    [t1,t2,m1,m2,~,~,~,~,~,~,evar,evardist] = estimate_vslite_params_v2_3(...
        squeeze(t_cal(q,:,gyrs{q})),squeeze(p_cal(q,:,gyrs{q})),...
        p_lat(q),TRW(gyrs{q},q)','nsamp',5000,'substep',1,'fourbetparams',fp);
    
    T1(q)=t1; T2(q)=t2;
    M1(q)=m1; M2(q)=m2;
    eVar(q)=evar;
    eVarDist{q}=evardist;
    
end

%------------------------
% Estimate MSE on my own
%------------------------

% Run VS-Lite
vsl_trw=zeros(pNum,length(c_yrs));
for q=1:pNum
    [vsl_trw(q,:)] = VSLite_v2_3(calPeriod(1),calPeriod(2),p_lat(q),T1(q),T2(q),...
        M1(q),M2(q),squeeze(t_cal(q,:,:)),squeeze(p_cal(q,:,:)),...
        'lbparams',[0.76,0.01,0.093,4.886,5.80,1000,0.2,1]);
end

% Compute MSE and correlation on tree ring overlap
MSE=zeros(pNum,1);
vsl_cr=zeros(pNum,1);
eVarUV=zeros(pNum,1);
for i=1:pNum
    vsl1=vsl_trw(i,gyrs{i});trw1=TRW(gyrs{i},i);
    MSE(i)=mse(vsl1(:),trw1(:));
    [~,~,eVarUV0]=regress(trw1(:),[ones(size(vsl1(:))) vsl1(:)]);
    eVarUV(i)=mean(eVarUV0.^2);
    vsl_cr(i)=nancorr(vsl1(:),trw1(:));
end


end

