
% calculates a 2-way least squares fit from weighted data
% see http://eos.whoi.edu/12.747/notes/lect03/lectno03.html

close all
clear all

load actXLSData_5Days

cw=0.4; % water leaving radiance

% put all tss, poc, bac, and fluor from each site into one array
tss=[cbl_tss;hawaii_tss;maine_tss;michigan_tss;michigan2_tss;mlml_tss;skio_tss;usf_tss];
tss_sd=[cbl_tss_sd;hawaii_tss_sd;maine_tss_sd;michigan_tss_sd;michigan2_tss_sd;mlml_tss_sd;skio_tss_sd;usf_tss_sd];
POC=[cbl_poc;hawaii_poc;maine_poc;michigan_poc;michigan2_poc;mlml_poc;skio_poc;usf_poc];
POC_sd=[cbl_poc_sd;hawaii_poc_sd;maine_poc_sd;michigan_poc_sd;michigan2_poc_sd;mlml_poc_sd;skio_poc_sd;usf_poc_sd];
bac=[cbl_bac4h;hawaii_bac4h;maine_bac4h;michigan_bac4h;michigan2_bac4h;mlml_bac4h;skio_bac4h;usf_bac4h]-cw;
fluor=[cbl_fluor;hawaii_fluor;maine_fluor;michigan_fluor;michigan2_fluor;mlml_fluor;skio_fluor;usf_fluor];

% put all turbidity data from each site into arrays by instrument
wetlabs=[cbl_wetlabs;hawaii_wetlabs;maine_wetlabs;michigan_wetlabs;michigan2_wetlabs;mlml_wetlabs;skio_wetlabs;usf_wetlabs];
ysi=[cbl_ysi;hawaii_ysi;maine_ysi;michigan_ysi;michigan2_ysi;mlml_ysi;skio_ysi;usf_ysi];
insitu=[cbl_insitu;hawaii_insitu;maine_insitu;michigan_insitu;michigan2_insitu;mlml_insitu;skio_insitu;usf_insitu];
mcvan=[cbl_mcvan;hawaii_mcvan;maine_mcvan;michigan_mcvan;michigan2_mcvan;mlml_mcvan;skio_mcvan;usf_mcvan];
aqua=[cbl_aqua;hawaii_aqua;maine_aqua;michigan_aqua;michigan2_aqua;mlml_aqua;skio_aqua;usf_aqua];

% make an index array so we can remove sites in some analyses (like Hawaii)
% index 1 (records 1-50): cbl (1-15 for 5 days)
% index 2 (next 71 records): hawaii (16-29 for 5 days)
% index 3 (next 63 records): maine (30-41 for 5 days)
% index 4 (next 59 records): michigan (41-54 for 5 days)
% index 5 (next 52 records): michigan2 (54-67 for 5 days)
% index 6 (next 80 records): mlml (67-82 for 5 days)
% index 7 (next 49 records): skio (82-91 for 5 days)
% index 8 (next 39 records): usf (91-91 for 5 days)

% index=[ones(50,1)*1;ones(71,1)*2;ones(63,1)*3;ones(59,1)*4;ones(52,1)*5;ones(80,1)*6;ones(49,1)*7;ones(39,1)*8]; % for all data, not just first 5 days
index=[ones(15,1)*1;ones(13,1)*2;ones(11,1)*3;ones(13,1)*4;ones(13,1)*5;ones(15,1)*6;ones(9,1)*7;ones(1,1)*8]; % first 5 days

% remove Hawaii
I=find(index~=2);

% X=tss or POC
% XP=std in tss or POC (tss_sd or POC_sd)
% Y=optical variable (mcvan, wetlabs)
% YP ones(1:length(optical data))*std of instrument e.g. 0.04*signal or a constant
  % signal noise, SDVB and Formazin added:
  % wetlabs: 100+(4200/85)*Y and 100+(4200/85)*Y;
  % aquatec: 0.02+(0.15/400)*Y and 0.02+(1.3/260)*Y;
  % insitu: 0.02+(0.24/580)*Y and 0.02+(0.08/270)*Y;
  % mcvan: 0.3+(0.55/83)*Y and 0.3+(0.80/100)*Y;
  % ysi: 0.02+(0.35/580)*Y and 0.05+(1.05/250)*Y;
% tl set as 10^-3.
% make sure to plot data and line, e.g.:
% plot(x,y,'*',m*x+b,'-')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO STEPS 1, 2, 3, AND 4
%%%%%%%%%%%%%%%%%%%%%%%% (STEP 1) select X, XP for a site:
% X=tss;
% XP=tss_sd;

% X=tss(I);
% XP=tss_sd(I);

% X=cbl_tss;
% XP=cbl_tss_sd;
% XP=min(cbl_tss_sd)*ones(1,length(cbl_tss))';

% X=hawaii_tss;
% XP=hawaii_tss_sd;

% X=maine_tss;
% XP=maine_tss_sd;

% X=michigan_tss;
% XP=michigan_tss_sd;

% X=michigan2_tss;
% XP=michigan2_tss_sd;

% X=mlml_tss;
% XP=mlml_tss_sd;

% X=skio_tss;
% XP=skio_tss_sd;

% X=usf_tss;
% XP=usf_tss_sd;

%%%%%%%%%%%%%%%%%%%%%%%% (STEP 2) select instrument by site, one at a time:
% Y=wetlabs;
% Y=aqua;
% Y=insitu;
% Y=mcvan;
% Y=ysi;
% Y=bac;

% Y=wetlabs(I);
% Y=aqua(I);
% Y=insitu(I);
% Y=mcvan(I);
% Y=ysi(I);
% Y=bac(I);

% Y=cbl_wetlabs;
% Y=cbl_aqua;
% Y=cbl_insitu;
% Y=cbl_mcvan;
% Y=cbl_ysi;
% Y=cbl_bac;

% Y=hawaii_wetlabs;
% Y=hawaii_aqua;
% Y=hawaii_insitu;
% Y=hawaii_mcvan;
% Y=hawaii_ysi;
% Y=hawaii_bac;

% Y=maine_wetlabs;
% Y=maine_aqua;
% Y=maine_insitu;
% Y=maine_mcvan;
% Y=maine_ysi;
% Y=maine_bac;

% Y=michigan_wetlabs;
% Y=michigan_aqua;
% Y=michigan_insitu;
% Y=michigan_mcvan;
% Y=michigan_ysi;
% Y=michigan_bac;

% Y=michigan2_wetlabs;
% Y=michigan2_aqua;
% Y=michigan2_insitu;
% Y=michigan2_mcvan;
% Y=michigan2_ysi;
% Y=michigan2_bac;

% Y=mlml_wetlabs;
% Y=mlml_aqua;
% Y=mlml_insitu;
% Y=mlml_mcvan;
% Y=mlml_ysi;
% Y=mlml_bac;

% Y=skio_wetlabs;
% Y=skio_aqua;
% Y=skio_insitu;
% Y=skio_mcvan;
% Y=skio_ysi;
% Y=skio_bac;

% Y=usf_wetlabs;
% Y=usf_aqua;
% Y=usf_insitu;
% Y=usf_mcvan;
% Y=usf_ysi;
% Y=usf_bac;

%%%%%%%%%%%%%%%%%%%%%%%% (STEP 3) select instrument error:
% YP=3*80/4000*ones(1,length(Y))'; % wetlabs std
% YP=0.02+(0.15/400)*Y; % aquatec std
% YP=0.02+(0.24/580)*Y; % insitu std
% YP=0.3+(0.55/83)*Y; % mcvan std
% YP=0.02+(0.35/580)*Y; % ysi std
% YP=0.01+2/30*Y % bac std 

tl=10^-3; % adjust this value to -2, then -1 if line is bad
siteID=[1];

% remove NaNs
I1=find(isfinite(X.*Y));

if isempty(I1)==1
    disp('no data!!!')
    break
end

%[MC,b3,r3,sm3,sb3,xbar3,ybar3] = lsqfitma(X(I1),Y(I1));
%med_dev=median(abs(Y(I1)-(b3+MC*(X(I1)))));

[mz,bz,rz,smz,sbz,xz,yz] = lsqfityw(Y(I1),X(I1),XP(I1));
med_dev=median(abs(Y(I1)-(-bz/mz+1/mz*(X(I1)))));
    
%check line
figure(1); plot(X(I1),Y(I1),'*',X(I1),X(I1)/mz+-bz/mz,'-');
hold on
    for i=1:length(I1)
       plot([X(I1(i)),X(I1(i))],[Y(I1(i))-YP(I1(i)),Y(I1(i))+YP(I1(i))]) % error bar left
       plot([X(I1(i))-XP(I1(i)),X(I1(i))+XP(I1(i))],[Y(I1(i)),Y(I1(i))]) % error bar right
       plot([X(I1(i))-XP(I1(i)),X(I1(i))-XP(I1(i))],[Y(I1(i))-0.3,Y(I1(i))+0.3]) % end bracket
       plot([X(I1(i))+XP(I1(i)),X(I1(i))+XP(I1(i))],[Y(I1(i))-0.3,Y(I1(i))+0.3]) % end bracket
    end
I=find(isfinite(X.*Y) & abs(Y-(-bz/mz+1/mz*(X)))<3*abs(med_dev));
figure(1); plot(X(I),Y(I),'r*');
hold off

n_outliers_removed=length(I1)-length(I);
if length(I)~=length(I1)
    disp('filtered out outliers')
    disp(length(I1)-length(I))
end

n=length(I);
[m4,b4,r4,sm4,sb4] = lsqfitgm(X(I),Y(I)); %recalculating regression w/o outliers
figure(1); 
hold on;
plot(X(I1),Y(I1),'*',X(I1),X(I1)*m4+b4,'k-');
hold off

M2=[m4 b4 sm4 sb4 r4];
[m5,b5,r5,sm5,sb5,x5,y5] = lsqfityw(Y(I),X(I),XP(I));
M3=[1/m5 -b5/m5 sm5/m5 sb5/m5 r5];

% get tss correlation coefficient
a=corr(X(I),Y(I));
z=0.5*log((1+a)/(1-a));
v=1/(n-3);
alpha=exp(2*(z-sqrt(v)*1.96));
beta=exp(2*(z+sqrt(v)*1.96));
corr_min_tss=(alpha-1)/(alpha+1);
corr_max_tss=(beta-1)/(beta+1);
    
C=[a corr_min_tss corr_max_tss];

% get tss rank correlation coefficient
[r,t,p] = spear(X(I),Y(I));
z3=0.5*log((1+r)/(1-r));
v3=1/(n-3);
alpha3=exp(2*(z3-sqrt(v3)*1.96));
beta3=exp(2*(z3+sqrt(v3)*1.96));
rank_corr_min_tss=(alpha3-1)/(alpha3+1);
rank_corr_max_tss=(beta3-1)/(beta3+1);
    
RC=[r rank_corr_min_tss rank_corr_max_tss];

if isfinite(I) % if we have data
    xlswrite('twoWayLsf.xls', siteID, 'tss (all)', 'A2');
    xlswrite('twoWayLsf.xls', C, 'tss (all)', 'H2:J2');
    xlswrite('twoWayLsf.xls', RC, 'tss (all)', 'K2:M2');
    xlswrite('twoWayLsf.xls', n_outliers_removed, 'tss (all)', 'N2');
    xlswrite('twoWayLsf.xls', tl, 'tss (all)', 'O2');
    xlswrite('twoWayLsf.xls', M2, 'tss (all)', 'Q2:U2');
    xlswrite('twoWayLsf.xls', M3, 'tss (all)', 'V2:Z2');
    
    % call lsqcubic function to calculate a 2-way least squares fit from weighted data
    [m,b,rw,sm,sb,xbar,ybar,ct]=lsqcubic(X(I),Y(I),XP(I),YP(I),tl);
    M=[length(I) m b sm sb rw];
    
    % check plots
    figure(2); plot(X(I),Y(I),'*',X(I),m*X(I)+b,'-');
    hold on
    for i=1:length(I1)
       plot([X(I1(i)),X(I1(i))],[Y(I1(i))-YP(I1(i)),Y(I1(i))+YP(I1(i))])
       plot([X(I1(i))-XP(I1(i)),X(I1(i))+XP(I1(i))],[Y(I1(i)),Y(I1(i))])
       plot([X(I1(i))-XP(I1(i)),X(I1(i))-XP(I1(i))],[Y(I1(i))-0.3,Y(I1(i))+0.3])
       plot([X(I1(i))+XP(I1(i)),X(I1(i))+XP(I1(i))],[Y(I1(i))-0.3,Y(I1(i))+0.3])
    end
    hold off
    
    xlswrite('twoWayLsf.xls', ct, 'tss (all)', 'P2');
    xlswrite('twoWayLsf.xls', M, 'tss (all)', 'B2:G2');
else
    disp('no data!!!!')
end

%%%%%%%%%%%%%%%%%%%%%%%% (STEP 4) open twoWayLsf.xls and insert a row below header
%%%%%%% run again with new variables 
   