function [maxCos, cos, clusterID, totScore] = QAscores_landsat8_matrix(test_Rrs, test_lambda)
% Quality assurance system for Rrs spectra (version 1)
%
% Author: Jianwei Wei, University of Massachusetts Boston
% Email: Jianwei.Wei@umb.edu
% Nov-01-2016
% updated March-28-2017 for Landsat 8
% updated April-12-2017; following QAscores_matrix.m
%
% ------------------------------------------------------------------------------
% KNOWN VARIABLES :   ref_nRrs   -- Normalized Rrs spectra per-determined from water clustering (23x9 matrix)  
%                     ref_lambda -- Wavelengths for ref_nRrs (1x9 matrix)
%                     upB        -- Upper boundary (23x9 matrix)
%                     lowB       -- Lower boundary (23x9 matrix)
%
% INPUTS:            test_Rrs   -- Rrs spectra for testing (units: sr^-1); a row vector
%                    test_lambda-- Wavelengths for test_Rrs
%
% OUTPUTS:  maxCos     -- maximum cosine values
%           cos        -- cosine values for every ref_nRrs spectra
%           clusterID  -- idenfification of water types (from 1-23)
%           totScore   -- total score assigned to test_Rrs
% ------------------------------------------------------------------------------
% 
% NOTE:
%         1) Nine wavelengths [443, 482, 561, 655] are assumed in the model
%         2) If your Rrs data were measured at other wavelength, e.g. 440nm, you may want to change 440 to 443 before the model run;
%             or modify the code below to find a cloest wavelength from the nine bands.
%         3) The latest version may be found online at HTTP://oceanoptics.umb.edu/score_metric
%
% Reference:
%         1)  Wei, Jianwei; Lee, Zhongping; Shang, Shaoling (2016). A system
%             to measure the data quality of spectral remote sensing
%             reflectance of aquatic environments. Journal of Geophysical Research, 
%             121, doi:10.1002/2016JC012126
%         2)  Wei, Jianwei; Lee, Zhongping(2017). Cross-platform validation and assessment of
%             Landsat-8 ocean color measruements in coastal turbid waters and coral reefs. 
%             Remote Sensing of Environment,
%             ####, ####
%    

%% check input data
[row_lam, len] = size(test_lambda);
if( row_lam ~= 1 )
    test_lambda = test_lambda';
    [row_lam, len] = size(test_lambda);
end

[row, col] = size(test_Rrs);
if( len~=col && len~=row)
    error('Rrs and lambda size mismatch, please check the input data!');
elseif( len == row )
    test_Rrs = test_Rrs';
end

%% 
ref_lambda = [... 
4.4300000e+02   4.8200000e+02   5.6100000e+02   6.5500000e+02];
ref_nRrs = [...
   8.3833649e-01   5.3317452e-01   1.0410211e-01   1.3057839e-02
   7.9842099e-01   5.8424520e-01   1.4234334e-01   1.9153490e-02
   7.6095686e-01   6.2017453e-01   1.8621319e-01   2.7343167e-02
   7.0260970e-01   6.5317576e-01   2.7612731e-01   4.5705366e-02
   6.4969075e-01   6.9266426e-01   3.0638474e-01   4.8486000e-02
   6.4907198e-01   6.6016554e-01   3.6916493e-01   6.8897188e-02
   5.8905023e-01   6.7971948e-01   4.3006590e-01   7.0411806e-02
   5.6775833e-01   6.4773550e-01   4.9661546e-01   9.9217817e-02
   5.2171284e-01   6.6562539e-01   5.2386109e-01   9.1983127e-02
   4.9976668e-01   6.1643112e-01   5.9347791e-01   1.2860846e-01
   5.2011932e-01   5.7677225e-01   6.0001790e-01   1.8048702e-01
   4.4126564e-01   5.9900177e-01   6.5007060e-01   1.4407249e-01
   4.4126038e-01   5.5338691e-01   6.7110426e-01   2.1325622e-01
   3.8460577e-01   5.2816131e-01   7.0698153e-01   2.6595912e-01
   3.6121112e-01   5.3769709e-01   7.3717604e-01   1.8089482e-01
   4.0918461e-01   4.9812757e-01   7.2199852e-01   2.3296006e-01
   3.5814155e-01   4.7442951e-01   7.0391624e-01   3.8492051e-01
   3.1076205e-01   4.6734058e-01   7.6256239e-01   3.1727574e-01
   3.1283767e-01   4.0538904e-01   6.8352286e-01   5.1548610e-01
   2.6721915e-01   4.1914901e-01   8.2571824e-01   2.5411151e-01
   2.5085988e-01   3.7063636e-01   7.8518970e-01   4.2082213e-01
   1.9490723e-01   3.1731938e-01   7.5752491e-01   5.3272522e-01
   1.7827700e-01   2.5261352e-01   7.3077695e-01   5.9592089e-01];
upB = [...
   8.8077395e-01   5.7366649e-01   1.3979904e-01   7.2643928e-02
   8.2011231e-01   6.1101827e-01   1.7668095e-01   4.8964652e-02
   7.8312650e-01   6.4397993e-01   2.4554703e-01   1.1119774e-01
   7.3680439e-01   6.7505872e-01   3.3724260e-01   9.4524130e-02
   6.8488196e-01   7.2737735e-01   3.5140255e-01   8.5595775e-02
   6.7963122e-01   6.8983117e-01   4.3312213e-01   1.6358164e-01
   6.1777156e-01   7.0349931e-01   4.6951949e-01   1.0459642e-01
   6.0605250e-01   6.7630741e-01   5.5460376e-01   1.5681959e-01
   5.6585527e-01   7.0161677e-01   5.6711185e-01   1.3047263e-01
   5.4101332e-01   6.4753698e-01   6.5207457e-01   1.8861404e-01
   5.6441964e-01   6.1848847e-01   6.7008393e-01   2.6657980e-01
   4.7768742e-01   6.8667756e-01   7.0148648e-01   2.0794128e-01
   4.7693577e-01   5.8116510e-01   7.4229731e-01   2.9279236e-01
   4.2351320e-01   5.6498537e-01   7.5404612e-01   3.3064146e-01
   4.0343328e-01   5.9950669e-01   7.8554540e-01   2.5276890e-01
   4.5577661e-01   5.4174497e-01   7.7807259e-01   4.7055969e-01
   4.0867966e-01   5.2049167e-01   7.4889201e-01   4.5670457e-01
   3.6769475e-01   5.3337054e-01   8.0371765e-01   3.9917274e-01
   3.6297084e-01   4.5859195e-01   7.4908314e-01   6.1114082e-01
   3.3444764e-01   4.9251118e-01   8.8552657e-01   3.7029967e-01
   3.1816068e-01   4.4742265e-01   8.4439398e-01   4.9922956e-01
   2.5306709e-01   3.6438024e-01   7.9089646e-01   5.9203932e-01
   2.3483514e-01   3.2898054e-01   8.2628841e-01   7.3360524e-01];
lowB = [...
   8.0877434e-01   4.6872142e-01   6.6966712e-02   4.4371662e-03
   7.7322322e-01   5.5976302e-01   1.1850104e-01   6.2171308e-03
   7.3579512e-01   5.9651685e-01   1.5957893e-01   1.2598273e-02
   6.7451468e-01   6.2673417e-01   2.2497670e-01   1.9270390e-02
   5.9686781e-01   6.7473551e-01   2.5550806e-01   2.7409519e-02
   6.1965571e-01   6.2032844e-01   3.2256814e-01   3.6681571e-02
   5.5003558e-01   6.4787983e-01   3.6967908e-01   4.5524034e-02
   5.3127430e-01   6.0126790e-01   4.6086274e-01   6.3882415e-02
   4.4645052e-01   6.3645307e-01   4.7093160e-01   3.2847153e-02
   4.6023688e-01   5.7629078e-01   5.5632897e-01   7.6892553e-02
   4.8538188e-01   5.2226266e-01   5.1339565e-01   1.3182796e-01
   4.0209761e-01   5.4999306e-01   5.6789526e-01   4.2944367e-02
   3.9010381e-01   5.0388361e-01   6.2154818e-01   1.3465256e-01
   3.3667364e-01   4.8734382e-01   6.6727370e-01   1.9689628e-01
   2.8143408e-01   4.9852582e-01   6.8435589e-01   6.2352051e-02
   3.3000436e-01   3.6595423e-01   6.8028500e-01   1.8473484e-01
   2.9516573e-01   4.3104477e-01   6.6575523e-01   3.3282180e-01
   2.5177392e-01   4.1654010e-01   7.2469690e-01   2.3479958e-01
   2.7315562e-01   3.3894210e-01   6.2587051e-01   4.6776411e-01
   1.9193855e-01   3.4983146e-01   7.8658023e-01   1.1111042e-01
   1.4348401e-01   2.8436176e-01   7.2388575e-01   3.0357637e-01
   1.4551303e-01   2.4965952e-01   7.3278833e-01   4.7533842e-01
   1.1234819e-01   1.7827135e-01   6.2463154e-01   4.6864217e-01];    
   
[refRow,refCol]=size(ref_nRrs);

%% match the ref_lambda and test_lambda
idx0 = []; % for ref_lambda 
idx1 = []; % for test_lambda

for i = 1 : length(test_lambda)
    pos = find(ref_lambda==test_lambda(i));
    if isempty(pos)
        idx1(i) = NaN;
    else
        idx0(i) = pos;
        idx1(i) = i;
    end
end

pos = isnan(idx1);  idx1(pos) = [];

test_lambda = test_lambda(idx1); test_Rrs = test_Rrs(:,idx1);
ref_lambda = ref_lambda(idx0); ref_nRrs = ref_nRrs(:,idx0); 
upB = upB(:,idx0); lowB = lowB(:,idx0); 

%% match the ref_nRrs and test_Rrs
% keep the original value
test_Rrs_orig = test_Rrs;

% inPosNan = isnan(test_Rrs);    %%%% Jianwei commented out these two lines
% test_Rrs(inPosNan) = 0;

% test_lambda(pos) = []; test_Rrs(:,pos) = [];
% ref_lambda(pos) = []; ref_nRrs(:,pos) = []; upB(:,pos) = []; lowB(:,pos) = []; 
   
%% nromalization
[inRow, inCol] = size(test_Rrs);

% transform spectrum to column, inCol*inRow
test_Rrs = test_Rrs';
test_Rrs_orig = test_Rrs_orig';

% inCol*inRow
nRrs_denom=sqrt(nansum(test_Rrs.^2));
nRrs_denom = repmat(nRrs_denom,[inCol,1]);
nRrs = test_Rrs./nRrs_denom;  

% SAM input, inCol*inRow*refRow 
test_Rrs2 = repmat(test_Rrs_orig,[1,1,refRow]);

%for ref Rrs, inCol*refRow*inRow 
test_Rrs2p = permute(test_Rrs2,[1,3,2]);

% pos = isnan(test_Rrs2);   %%%% Jianwei commented out these two lines
% test_Rrs2(pos) = 0;

% inCol*inRow*refRow  
nRrs2_denom=sqrt(nansum(test_Rrs2.^2));
nRrs2_denom = repmat(nRrs2_denom,[inCol,1]);
nRrs2 = test_Rrs2./nRrs2_denom;  
% inCol*refRow*inRow  
nRrs2 = permute(nRrs2,[1,3,2]);

%% adjust the ref_nRrs, according to the matched wavebands
%[row, ~] = size(ref_nRrs);

%%%% re-normalize the ref_adjusted
% for i = 1 : row
%     ref_nRrs_corr(i,:) = ref_nRrs(i,:)./sqrt(nansum(ref_nRrs(i,:).^2));  
% end
%
% transform spectrum to column, inCol*refRow
ref_nRrs = ref_nRrs';

% inCol*refRow*inRow 
ref_nRrs2 = repmat(ref_nRrs,[1,1,inRow]);

% % inCol*refRow*inRow                  %%%%%  Jianwei commented out
% pos = isnan(test_Rrs2p); 
% ref_nRrs2(pos) = 0;

% inCol*refRow*inRow 
ref_nRrs2_denom=sqrt(nansum(ref_nRrs2.^2));
ref_nRrs2_denom = repmat(ref_nRrs2_denom,[inCol,1]);
ref_nRrs_corr2 = ref_nRrs2./ref_nRrs2_denom;

%% Classification 
%%%% calculate the Spectral angle mapper
% for i = 1 : row    
%     cos(i,1) = nansum(ref_nRrs_corr(i,:).*nRrs)/...
%                 sqrt(nansum((ref_nRrs_corr(i,:)).^2) *nansum(nRrs.^2));    
% end
%cos = (ref_nRrs_corr'*nRrs)./sqrt(nansum(ref_nRrs_corr.^2)'*nansum(nRrs.^2)); 

% inCol*refRow*inRow 
cos_denom=sqrt(nansum(ref_nRrs_corr2.^2).*nansum(nRrs2.^2));
cos_denom = repmat(cos_denom,[inCol,1]);
cos = (ref_nRrs_corr2.*nRrs2)./cos_denom; 
% 1*refRow*inRow 
cos = sum(cos);
% refRow*inRow
cos = permute(cos,[2,3,1]);

% 1*inRow
[maxCos,clusterID] = max(cos);
posClusterID = isnan(maxCos);
%potential bug for vectorized code
%clusterID(pos) = NaN;

% if isnan(cos)
%     clusterID = NaN;
% else
%     clusterID = find(cos==maxCos);
% end

%% scoring
%C = [];
%[row, col] = size(ref_nRrs);

%%% re-normalzation
% for i = 1 : row
%     upB_corr(i,:) = upB(i,:)./sqrt(nansum(ref_nRrs(i,:).^2)); 
%     lowB_corr(i,:) = lowB(i,:)./sqrt(nansum(ref_nRrs(i,:).^2)); 
% end
% upB_corr = upB'./sqrt(nansum(ref_nRrs.^2)); 
% lowB_corr = lowB'./sqrt(nansum(ref_nRrs.^2)); 
upB_corr = upB'; 
lowB_corr = lowB'; 

%% comparison
% inCol*inRow
upB_corr2 = upB_corr(:,clusterID).*(1+0.005);
lowB_corr2 = lowB_corr(:,clusterID).*(1-0.005);
ref_nRrs2 = ref_nRrs(:,clusterID);

%normalization, take account of NaN bands
% ref_nRrs2(inPosNan') = 0;                     %%%%% Jianwei commented out
ref_nRrs2_denom=sqrt(nansum(ref_nRrs2.^2));
ref_nRrs2_denom = repmat(ref_nRrs2_denom,[inCol,1]);
upB_corr2 = upB_corr2 ./ ref_nRrs2_denom;
lowB_corr2 = lowB_corr2 ./ ref_nRrs2_denom;

upB_diff = upB_corr2 - nRrs;
lowB_diff = nRrs - lowB_corr2;

C = zeros(inCol,inRow);
pos = find( upB_diff>=0 & lowB_diff>=0 );
C(pos) = 1;

%process all NaN spectral 
C(:,posClusterID)=NaN;     %%%%%%% jianwei changed "0" to "NaN"

% for i = 1 : col
%     if nRrs(i)> upB_corr(clusterID,i)*(1+0.005) || ...
%                 nRrs(i)< lowB_corr(clusterID,i)*(1-0.005)
%         C(i,1) = 0;
%     else
%         C(i,1) = 1;
%     end 
% end

totScore = nanmean(C) ;