function [out] = outmeth1X(X, wave, tune, tuneR, xval, srdcv, sig, dispTF, srdoutTF, merit)
%This function does all merits for outlier detection and puts them into
%seperate blocks 1 for vector mean, 1 for each to space merit
%{
NOTE: Has all the values minmized (i.e. 1-cos, 1-r1)
    Inputs:
    X: spectra of data than is n samples by m wavelengths
    wave: the actual wavelengths measured in nm
		if unknown put 0
    tune: the tuning paramet you want to stop at for spatial merits
		if you want all put nan
	xval: a row spectrum to test if it is in calibration space
		if you want to do detection on the cal set put []
	srdcv: whether to do or not do internal srd cross-val
	dispTF: whether to display which merit is is calculating
	srdoutTF: whether to compute srdout


    Outputs:
    srdinsc: a cell array of merits, with the columns in each cell being
    sample out, and the rows being the merit (some corresponding to the number of
    eigenvectors used to build the space ), rows scalled by normallizing unit
    length

    srdin= a cell array of merits, with the columns in each cell being the
    sample out, and the rows being the merit (some corresponding to the number of
    eigenvectors used to build the space ) RAW VALUES

	NOTE: When you run a val. sample, my method still does a LOOCV on the
	calibration set and then at the end compares the validation sample to
	all the calibration samples
%}


scale='nul';


nMerits = length(merit);
srdin=cell(nMerits,1);
fusMerit=cell(nMerits,1);
fusMerit{1,1} = cell(1,1);
if nargin<9, srdoutTF=true; end
if nargin<8, dispTF=true; end
if nargin<7, sig=3; end
if nargin<6, srdcv=false; end
if nargin<5, xval=[]; end
if nargin<4, tuneR=[]; end
if nargin<3, tune=nan; end
if nargin<2, wave=0; end

out.input.X=X;
out.input.wave=wave;
out.input.tune=tune;
out.input.tuneR=tuneR;
out.input.xval=xval;
out.input.sig=sig;

i = 0;
j = 1;
%merit calculations
if sum(strcmpi(merit, 'vm'))
    i = i+1;
    if dispTF, disp('cos'), end
    srdin{i,1}(1,:)=1-bbcos(X, xval);
    fusMerit{1,1}{j,1}='cos'; j=j+1;
    
    if dispTF, disp('ed'), end
    srdin{i,1}(2,:)=bbed(X, xval);
    fusMerit{1,1}{j,1}='d_i'; j=j+1;
    
    if dispTF, disp('drt'), end
    srdin{i,1}(3,:)=bbdrt(X, xval);
     fusMerit{1,1}{j,1}='Det'; j=j+1;
     
    % Procrustus
    if dispTF, disp('PA'), end
    PAout=bbPA(X, xval);
    srdin{i,1}(4,:)=PAout.Ftm;
     fusMerit{1,1}{j,1}='T_i'; j=j+1;
     
    srdin{i,1}(5,:)=PAout.Htm;
    fusMerit{1,1}{j,1}='H_i'; j=j+1;
    
    srdin{i,1}(6,:)=PAout.ptm;
    fusMerit{1,1}{j,1}='p_i'; j=j+1;
    %making r1 and into vector to vector merits
    if dispTF, disp('r1'), end
    srdin{i,1}(7,:)=1-bbR1vv(X, xval);
    fusMerit{1,1}{j,1}='1-r_i'; j=j+1;
    if (wave~=0)
        %EISC/EISCd
        if dispTF, disp('EISCd'), end
        EISCDout=bbEISCD(X, wave, xval);
        srdin{1,1}(8,:)=EISCDout.nbtmd;
        fusMerit{1,1}{j,1}='EISCDb_1'; j=j+1;
        
        srdin{1,1}(9,:)=EISCDout.nXbtmd;
        fusMerit{1,1}{j,1}='EISCDXb_1'; j=j+1;
        
        srdin{1,1}(10,:)=EISCDout.nbtsod;
        fusMerit{1,1}{j,1}='EISCDb_2'; j=j+1;
        
        srdin{1,1}(11,:)=EISCDout.nXbtsod;
        fusMerit{1,1}{j,1}='EISCDXb_2'; j=j+1;
    end
end

k = 2;
%vector to space
if sum(strcmpi(merit, 'md'))
	i = i+1;
if dispTF, disp('MD'), end
srdin{i,1}=bbMD(X, tune, xval);
 fusMerit{k,1}='MD'; k=k+1;
end
if sum(strcmpi(merit, 'sint'))+sum(strcmpi(merit,'qres')) > 0
	i = i+1;
	if sum(strcmpi(merit, 'sint'))+sum(strcmpi(merit,'qres')) == 2
		i = i+1;
		if dispTF, disp('Sint/Qres'), end
		[srdin{i-1,1},srdin{i,1}]=bbSin(X, tune, xval);% sint in 3 (i-1), qres in 4 (i)
         fusMerit{k,1}='Sin'; k=k+1;
         fusMerit{k,1}='Q'; k=k+1;
	elseif sum(strcmpi(merit, 'sint')) 
		[srdin{i,1}, ~] = bbSin(X, tune, xval);
        fusMerit{k,1}='Sin'; k=k+1;
	elseif sum(strcmpi(merit, 'qres'))
		[~,srdin{i,1}] = bbSin(X, tune, xval);
        fusMerit{k,1}='Q'; k=k+1;
	end
end

if sum(strcmpi(merit, 'dc'))
	i = 1+i;
if dispTF, disp('dc'), end
srdin{i,1}=bbDC(X, tune, xval);
fusMerit{k,1}='DC';
end


% Remove rows in each block that have NaN's (occurs at high tune parms)
for i=1:nMerits
	A=srdin{i,1};
	A(any(isnan(A),2),:) = [];
	srdin{i,1}=A;
end

if ~isempty(tuneR) && sum(isnan(tuneR))>1
	srdin = inred(srdin, tuneR);
end

% scaling
srdinsc= scalingvs(srdin, scale, nMerits);

if srdoutTF
	srdout=srdPack(srdinsc, 'max');
	srdout.thres=thresCal(srdout, sig);
	
	if srdcv
		srdoutcv=srdPack(srdinsc, 'max', 'both', 'lmo', 50, 15, true);
	end
end
out.srdin=srdinsc;
out.srdinunsc=srdin;
if srdoutTF
	out.srdout=srdout;
	if srdcv
		out.srdoutcv=srdoutcv;
	end
end
out.fusMerit=fusMerit;

end

function outcos=bbcos(X, xval)
% Finds the cos of theta between a vector and the mean of a subset the rest of the
% spectra in X using cv
%{
        Inputs:
            X:  observed data that is n x m (samples x wavelengths)
            nIter: number of cross val iterations
            percent: a percentage going into the 'cal' set per iteration
        Outputs:
            coscv: vector (1,n) of cosine of angle between the sample and the mean of the
                        rest of X
        
%}
n=size(X,1);
if isempty(xval),
	outcos=nan(1,n);
else outcos=nan(1,n+1);
end
for i=1:n
	xso=X(i, :);
	Cal=X;
	Cal(i,:)=[];
	cal=mean(Cal,1);
	dot1=dot(xso, cal);
	normxso=norm(xso, 2);
	normcal=norm(cal,2);
	outcos(i)=dot1/normxso/normcal;
end % i loop
if ~isempty(xval)
	Cal=X;
	cal=mean(Cal,1);
	dot1=dot(xval,cal);
	normxval=norm(xval,2);
	normcal=norm(cal,2);
	outcos(n+1)=dot1/normxval/normcal;
	
end % if
end

function [ed] = bbed(X, xval)
%Finds the eucdean distanace between a vector and the mean of the rest of the spec in X
% Comparing each spectra to the rest of the spectra in X
%{
        Inputs:
            X:  observed data that is n x m (samples x wavelengths)
            xval: a validation sample (row vector)
        Outputs:
            edist: vector (1,n) of ed of comparing each specta to the rest of X
%}
n=size(X,1);
if isempty(xval),ed=nan(1,n);else ed=nan(1,n+1);end
for i=1:n
	xso=X(i,:)'; % column vector
	Cal=X;
	Cal(i,:)=[];
	cal=mean(Cal,1)'; % column vector
	ed(1,i)= norm(xso-cal,2);
end
if ~isempty(xval)
	Cal=X;
	cal=mean(Cal,1)';
	ed(1,n+1)=norm(xval'-cal,2);
end
end

function [drt] = bbdrt(X, xval)
%Finds the eucdean distanace between a vector and the mean of the rest of the spec in X
% Comparing each spectra to the rest of the spectra in X
%{
        Inputs:
            X:  observed data that is n x m (samples x wavelengths)
            xval: a validation sample (row vector)
        Outputs:
            edist: vector (1,n) of ed of comparing each specta to the rest of X
%}
n=size(X,1);
if isempty(xval),drt=nan(1,n);
else drt=nan(1,n+1);
end

for i=1:n
	xso=X(i,:); % row vector
	Cal=X;
	Cal(i,:)=[];
	cal=mean(Cal,1); % row vector
	Xprime=[cal; xso]; % augment sample out vector and vector mean
	drt(1,i)=det(Xprime*Xprime'); % a 2x2 matrix, if you do covariance (mxm) the drt is zero
end
if ~isempty(xval)
	Cal=X;
	cal=mean(Cal,1); % rowvector
	Xprime=[cal; xso];
	drt(1,n+1)=det(Xprime*Xprime');
	
end
end

function [out] = bbPA(X, xval)
%This function does procrustes anaylsis of X with unconstrained and
%contrained methods. This method uses rotation and streching to go from one
%'space' to another 'space' (with vm, the spaces are just 2-D)
%{
Note: This functions sets the target as both the sample out and as the mean
Note: This function used covariance arrays for the spectra (so and mean)
    Inputs:
        X:  observed data that is n x m (samples x wavelengths)
    Outputs:
        F: uncontrained transformation (1, n)
        H: contrained rotatoin transformation, shape measure (1,n)
        p: contrained dilation (same amount of streching throughout) transformation, magnitude measure (1,n)
%}
n=size(X,1);
% Preallocation
if isempty(xval),
	Ftm=zeros(1,n);
	Htm=zeros(1,n);
	ptm=zeros(1,n);
else
	Ftm=zeros(1,n+1);
	Htm=zeros(1,n+1);
	ptm=zeros(1,n+1);
end

for i=1:n
	xso=X(i,:); % row vector
	Cal=X;
	Cal(i,:)=[];
	
	cal=mean(Cal,1); % row vector
	covcal=(cal'*cal); %covaranince space
	covxso=(xso'*xso);
	
	%% Target is mean spectrum
	% Uncontrained PA
% 	[Um, Sm, Vm]= svd(covcal); %svd of covarriance matrix of cal bar (Y)
% 	[Uxso, Sxso, Vxso]= svd(covxso); %svd of covarriance matrix of  so (X)
    Um = cal./norm(cal);
    Vm = cal./norm(cal);
    Sm = norm(cal)^2;
    
    Uxso = xso./norm(xso);
    Vxso = xso./norm(xso);
    Sxso = norm(xso)^2;
	Ftt=(Vm(:,1)/(Sm(1,1))*Um(:,1)')*(covcal); % mean to mean  (F1)
	Fvt=(Vxso(:,1)/(Sxso(1,1))*Uxso(:,1)')*(covcal); % sample out to mean (F2)
	Ftm(1,i)=norm(Fvt-Ftt, 'fro'); % Frob. norm of F2-F1
	
	%% Contrained PA
	[Uh,Sh,Vh]=svds(covxso'*covca,1); % svd of XsoXbar for H2
	Hvt=Uh(:,1:1)*Vh(:,1:1)'; % H2
	Httold=Uxso(:,1)*Vxso(:,1)'; % H1 but incorrect
	Htt=Um(:,1)*Vm(:,1)' ;% H1 fixed 2/6/2015
	Htm(1,i)= norm(Hvt-Htt, 'fro'); % Frob norm of H2-H1
	ptm(1,i)=abs(trace (Sh)/trace(covxso*covxso')-trace(Sm)/trace(covcal*covcal'));% |p2-p1|
	
end
%% check validation sample
if ~isempty(xval)
	Cal=X;
	
	cal=mean(Cal,1);
	covcal=(cal'*cal);
	covval=(xval'*xval);
	
	%% Target is mean spectrum
	% Uncontrained PA
	[Um, Sm, Vm]= svd(covcal); %svd of covarriance matrix of (Y)
	[Uv, Sv, Vv]= svd(covval); %svd of covarriance matrix of (X)
	Ftt=(Vm(:,1)/(Sm(1,1))*Um(:,1)')*(covcal); % mean to mean (F1)
	Fvt=(Vv(:,1)/(Sv(1,1))*Uv(:,1)')*(covcal); % sample to mean (F2)
	Ftm(1,n+1)=norm(Fvt-Ftt, 'fro');
	
	%% Contrained PA
	[Uh,Sh,Vh]=svd(covxso'*covcal);
	Hvt=Uh(:,1:1)*Vh(:,1:1)'; % H2
	Htt=Uv(:,1)*Vv(:,1)'; % H1
	Htm(1,n+1)= norm(Hvt-Htt, 'fro'); % H2-H1
	ptm(1,n+1)=abs(trace (Sh)/trace(covval*covval')-trace(Sm)/trace(covcal*covcal'));
end
%% Outs
out.Ftm=Ftm;
out.Htm=Htm;
out.ptm=ptm;
out.Httold=Httold;
out.Htt=Htt;
end

function [out] = bbEISCD(X, wave, xval)
%This functions performs an EISC and EISC-d analysis of X, and outputs the
%merits discussed on 4/14
%{
Inputs
	X:  observed data that is n x p (samples x wavelengths)
	wave: vector of wavelengths
	xval: row vector to do quality control
%}
[n,p]=size(X);
if iscolumn(wave), wave=wave';end
if isempty(xval),
	nbtmd=nan(1,n);
	nXbtmd=nan(1,n);
	nbtsod=nan(1,n);
	nXbtsod=nan(1,n);
else
	nbtmd=nan(1,n+1);
	nXbtmd=nan(1,n+1);
	nbtsod=nan(1,n+1);
	nXbtsod=nan(1,n+1);
end


wave=wave';
wave=wave./1000;
wave2=wave.^2;
lnwave=log(wave); %natural log
one=ones(p,1);

for i=1:n % num of samples in X
	xso=X(i,:)'; % column vector
	Cal=X;
	Cal(i,:)=[];
	cal=mean(Cal,1)'; % column vector
	
	% EISC-d
	% target= mean spectrum (cal)
	xt=cal;
	xl=xso;
	x2=xl.^2;
	der1=diff(xl);
	der2=diff(xl,2);
	z=p-1;
	y=xt-xl;
	
	X2=[one(2:z,1),xl(2:z,1),x2(2:z,1),der1(2:z,1),der2,wave(2:z,1),wave2(2:z,1),lnwave(2:z,1)];
	btmd= (X2'*X2)\X2'*y(2:end-1,:);
	nbtmd(1,i)=norm(btmd);
	Xb=X2*btmd;
	nXbtmd(1,i)=norm(Xb,2);
	
	% target= mean spectrum
	xt=xso;
	xl=cal;
	x2=xl.^2;
	der1=diff(xl);
	der2=diff(xl,2);
	z=p-1;
	y=xt-xl;
	
	X2=[one(2:z,1),xl(2:z,1),x2(2:z,1),der1(2:z,1),der2,wave(2:z,1),wave2(2:z,1),lnwave(2:z,1)];
	btsod= (X2'*X2)\X2'*y(2:end-1,:);
	nbtsod(1,i)= norm(btsod,2);
	Xb=X2*btsod;
	nXbtsod(1,i)=norm(Xb,2);
end

%% for validation sample
if ~isempty(xval)
	xval=xval';% column vector
	Cal=X;
	cal=mean(Cal,1)'; % column vector
	
	%%EISC-d
	% target= mean spectrum (cal)
	xt=cal;
	xl=xval;
	x2=xl.^2;
	der1=diff(xl);
	der2=diff(xl,2);
	z=p-1;
	y=xt-xl;
	
	X2=[one(2:z,1),xl(2:z,1),x2(2:z,1),der1(2:z,1),der2,wave(2:z,1),wave2(2:z,1),lnwave(2:z,1)];
	btmd= (X2'*X2)\X2'*y(2:end-1,:);
	nbtmd(1, n+1)=norm(btmd);
	Xb=X2*btmd;
	nXbtmd(1, n+1)=norm(Xb,2);
	
	% target= sammple out(val)
	xt=xval;
	xl=cal;
	x2=xl.^2;
	der1=diff(xl);
	der2=diff(xl,2);
	z=p-1;
	y=xt-xl;
	
	X2=[one(2:z,1),xl(2:z,1),x2(2:z,1),der1(2:z,1),der2,wave(2:z,1),wave2(2:z,1),lnwave(2:z,1)];
	btsod= (X2'*X2)\X2'*y(2:end-1,:);
	nbtsod(1, n+1)= norm(btsod,2);
	Xb=X2*btsod;
	nXbtsod(1, n+1)=norm(Xb,2);
end
%% Outputs


out.nbtmd=nbtmd;
out.nXbtmd=nXbtmd;

out.nbtsod=nbtsod;
out.nXbtsod=nXbtsod;
end

function [md] = bbMD(X, tune, xval)
%This function does finds the Mahalanobis distance between a sample out and
%the centroid of a space.
%{
Inputs:
	X:  observed data that is n x m (samples x wavelengths)
	xval: sample in question (default [])
	tune: truncates additional eigs (high ones), if blank the default is the rank of
	Cal (defalut is nan)
Outputs:
	md: matrix of Mahalanobis distances that has each column as a
	sample out and each row has the number of eigenvectors used to span
	X'X
      
%}
n=size(X,1);
if isempty(xval),
	md=NaN(n-1,n);
else
	md=NaN(n,n+1);
end

for i=1:n
	xso=X(i,:)'; % column vector, sample out
	Cal=X;
	Cal(i,:)=[];
	calmean=mean(Cal,1)'; % column mean vector of library
	[Calmc, ~]=meancent(Cal);
	covCal=Calmc'*Calmc; %covariance matrix of cal space
	covCal=covCal./(size(Cal,1)); % need to scale so covariance is the same magnitude
	rk=rank(covCal);
	[U,S,V]=svd(covCal);
	
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for j=1:tune
		calcovinv=V(:, 1:j)/(S(1:j,1:j))*U(:,1:j)';
		md(j,i)=sqrt((xso-calmean)'*calcovinv*(xso-calmean));
	end
end
%% for validation sample
if ~isempty(xval)
	Cal=X;
	calmean=mean(Cal,1)'; % column mean vector of library
	[Calmc, ~]=meancent(Cal);
	covCal=Calmc'*Calmc; %covariance matrix of cal space
	covCal=covCal./(size(Cal,1)); % need to scale so covariance is the same magnitude
	[U,S,V]=svd(covCal);
	
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for j=1:tune
		calcovinv=V(:, 1:j)/(S(1:j,1:j))*U(:,1:j)';
		md(j,n+1)=sqrt((xval'-calmean)'*calcovinv*(xval'-calmean));
	end
	
end
end % function

function [outSin, Qres] = bbSin(X, tune, xval)
%This function does finds the sine of theta a sample out and
%a space by using the projections, as well as the x residual (x projection)
%{
    Inputs:
        X:  observed data that is n x m (samples x wavelengths)
    Outputs:
        sin: matrix of sin theatas that has each column as a
        sample out and each row has the number of eigenvectors used to span
        Xcal
        Qres: the squared 2 norm of the orthogonal projection of each
        sample, equal to sqrt(xres)
%}
[n, p]=size(X);
if isempty(xval)
	outSin=NaN(n,n);Qres=NaN(n,n);
else
	outSin=NaN(n,n+1);Qres=NaN(n,n+1);
end
for i=1:n
	xso=X(i,:)'; %column vector, sample out
	Cal=X;
	Cal(i,:)=[]; % makes library matrix
	
	rk=rank(Cal);
	[~,~,V]=svd(Cal);
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for j=1:tune
		P=V(:,1:j)*V(:,1:j)';
		I=eye(p);
		valstar=(I-P)*xso;
		outSin(j,i)=norm(valstar,2)/norm(xso,2);
		Qres(j,i)=(outSin(j,i)*norm(xso,2))^2;
	end
	
end
%% for validation sample
if ~isempty(xval)
	Cal=X;
	
	[~,~,V]=svd(Cal);
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for j=1:tune
		P=V(:,1:j)*V(:,1:j)';
		I=eye(p);
		valstar=(I-P)*xval';
		outSin(j, n+1)=norm(valstar,2)/norm(xval,2);
		Qres(j,n+1)=(outSin(j,n+1)*norm(xval,2))^2;
	end
end
end % function

function [r1] = bbR1(X, tune, xval)
%This function finds the inner product colrreation  (r1) between two covariance
%spaces, one being the sample out and the other being the Cal set

[n, ~]=size(X);
if isempty(xval)
	r1=NaN(n,n);
else
	r1=NaN(n,n+1);
end

for i=1:n
	xso=X(i,:)'; % column vector
	Cal=X;
	Cal(i,:)=[];
	covxso=xso*xso';
	covCal=Cal'*Cal./(size(Cal,1));
	rk=rank(covCal);
	[Uso, Sso, Vso]=svd(covxso);
	[Uc, Sc, Vc]=svd(covCal);
	
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for j=1:tune
		num=trace(Vso(:,1)*(Sso(1,1))*Uso(:,1)'*Uc(:,1:j)*Sc(1:j,1:j)*Vc(:,1:j)');
		dem=sqrt(trace(Sso(1,1)^2)*trace(Sc(1:j,1:j)'*Sc(1:j, 1:j)));
		r1(j,i)=(num)/(dem);
	end % # of eigs
end % sample out

% for validation sample
if ~isempty(xval)
	Cal=X;
	
	covval=xval'*xval;
	covCal=Cal'*Cal./(size(Cal,1));
	[Uv, Sv, Vv]=svd(covval);
	[Uc, Sc, Vc]=svd(covCal);
	
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for k=1:tune
		num=trace(Vv(:,1)*(Sv(1,1))*Uv(:,1)'*Uc(:,1:k)*Sc(1:k,1:k)*Vc(:,1:k)');
		dem=sqrt(trace(Sv(1,1)^2)*trace(Sc(1:k,1:k)'*Sc(1:k,1:k)));
		r1(k, n+1)=(num)/(dem);
	end % # of eigs
end % if xval
end % function

function [r1] = bbR1vv(X, xval)
%This function finds the inner product colrreation  (r1) between two covariance
%spaces, one being the sample out and the other being the Cal set

[n, ~]=size(X);
if isempty(xval)
	r1=NaN(1,n);
else
	r1=NaN(1,n+1);
end

for i=1:n
	xso=X(i,:)'; % column vector
	Cal=X;
	Cal(i,:)=[];
	covxso=xso*xso';
	cal=mean(Cal,1)'; % column vector
	covcal=cal*cal';
	[Uso, Sso, Vso]=svd(covxso);
	[Uc, Sc, Vc]=svd(covcal);
	
	num=trace(Vso(:,1)*(Sso(1,1))*Uso(:,1)'*Uc(:,1:1)*Sc(1:1,1:1)*Vc(:,1:1)');
	dem=sqrt(trace(Sso(1,1)^2)*trace(Sc(1:1,1:1)'*Sc(1:1, 1:1)));
	r1(1,i)=(num)/(dem);
end % sample out

% for validation sample
if ~isempty(xval)
	Cal=X;
	cal=mean(Cal,1)'; % column vector
	covval=xval'*xval;
	covcal=cal*cal';
	[Uv, Sv, Vv]=svd(covval);
	[Uc, Sc, Vc]=svd(covcal);
	
	num=trace(Vv(:,1)*(Sv(1,1))*Uv(:,1)'*Uc(:,1:1)*Sc(1:1,1:1)*Vc(:,1:1)');
	dem=sqrt(trace(Sv(1,1)^2)*trace(Sc(1:1,1:1)'*Sc(1:1,1:1)));
	r1(1, n+1)=(num)/(dem);
	
end % if xval
end % function

function [DC] = bbDC(X, tune, xval)
%This functions finds the Divergence Criteria between 2 covariance spaces,
%one is the sample out, the other is the Cal Set (all the other samples but
%the sample out) TAKES ABSOULTE VALUE OF DC VALUES

n=size(X,1);
if isempty(xval);
	DC=NaN(n-1,n);
else
	DC=NaN(n,n+1);
end

for i=1:n
	xso=X(i,:)'; %column vector
	Cal=X;
	Cal(i,:)=[];
	cal=mean(Cal,1)'; % column vector
	covxso=xso*xso'; %S1
	covCal=Cal'*Cal./(size(Cal,1)); %S2
	rk=rank(covCal);
	[Uso, Sso, Vso]=svd(covxso);
	[Uc, Sc, Vc]=svd(covCal);
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	for j=1:tune
		num1=.5*trace((Uso(:,1)*Sso(1,1)*Vso(:,1)'-Uc(:,1:j)*Sc(1:j,1:j)*Vc(:,1:j)')*(Vc(:,1:j)/Sc(1:j,1:j)*Uc(:,1:j)'-Vso(:,1)/Sso(1,1)*Uso(:,1)'));
		%1/2*trace((S1-S2)(S2inv-S1inv))
		num2=.5*trace( (Vso(:,1)/Sso(1,1)*Uso(:,1)'+Vc(:,1:j)/Sc(1:j,1:j)*Uc(:,1:j)')*(xso-cal)*(xso-cal)');
		%1/2*(xso-xmean)'(S1inv+S2inv)(xso-xmean)
		DC(j,i)=abs(num1+num2);
	end % rk
	
end % n samples
% for validation sample
if ~isempty(xval)
	Cal=X;
	cal=mean(Cal,1)'; % column vector
	covval=xval'*xval; %S1
	covCal=Cal'*Cal./(size(Cal,1)); %S2
	[Uv, Sv, Vv]=svd(covval);
	[Uc, Sc, Vc]=svd(covCal);
	
	if isnan(tune)==true, tune=rk; end
	if tune>=rk, tune=rk; end
	
	for j=1:tune
		num1=.5*trace((Uv(:,1)*Sv(1,1)*Vv(:,1)'-Uc(:,1:j)*Sc(1:j,1:j)*Vc(:,1:j)')*...
			(Vc(:,1:j)/Sc(1:j,1:j)*Uc(:,1:j)'-Vv(:,1)/Sv(1,1)*Uv(:,1)'));
		%1/2*trace((S1-S2)(S2inv-S1inv))
		num2=.5*trace((Vv(:,1)/Sv(1,1)*Uv(:,1)'+Vc(:,1:j)/Sc(1:j,1:j)*Uc(:,1:j)')*(xval'-cal)*(xval'-cal)');
		%1/2*(S1inv+S2inv)(xso-xmean)'*(xso-xmean)
		DC(j,n+1)=abs(num1+num2);
	end % rk
end % if statement
end % function

function [DC] = bbDCvv(X, xval)
%This functions finds the Divergence Criteria between 2 covariance spaces,
%one is the sample out, the other is the Cal Set (all the other samples but
%the sample out) TAKES ABSOULTE VALUE OF DC VALUES

n=size(X,1);
if isempty(xval);
	DC=NaN(1,n);
else
	DC=NaN(1,n+1);
end

for i=1:n
	xso=X(i,:)'; %column vector
	Cal=X;
	Cal(i,:)=[];
	cal=mean(Cal,1)'; % column vector
	covxso=xso*xso'; %S1
	covcal=cal*cal'; %S2
	
	[Uso, Sso, Vso]=svd(covxso);
	[Uc, Sc, Vc]=svd(covcal);
	
	num1=.5*trace((Uso(:,1)*Sso(1,1)*Vso(:,1)'-Uc(:,1:1)*Sc(1:1,1:1)*Vc(:,1:1)')*...
		(Vc(:,1:1)/Sc(1:1,1:1)*Uc(:,1:1)'-Vso(:,1)/Sso(1,1)*Uso(:,1)'));
	%1/2*trace((S1-S2)(S2inv-S1inv))
	num2=.5*trace( (Vso(:,1)/Sso(1,1)*Uso(:,1)'+Vc(:,1:1)/Sc(1:1,1:1)*Uc(:,1:1)')*(xso-cal)*(xso-cal)');
	%1/2*(xso-xmean)'(S1inv+S2inv)(xso-xmean)
	DC(1,i)=abs(num1+num2);
	
end % n samples

% for validation sample
if ~isempty(xval)
	Cal=X;
	cal=mean(Cal,1)'; % column vector
	covval=xval'*xval; %S1
	covcal=cal*cal'; %S2
	[Uv, Sv, Vv]=svd(covval);
	[Uc, Sc, Vc]=svd(covcal);
	
	num1=.5*trace((Uv(:,1)*Sv(1,1)*Vv(:,1)'-Uc(:,1:1)*Sc(1:1,1:1)*Vc(:,1:1)')*...
		(Vc(:,1:1)/Sc(1:1,1:1)*Uc(:,1:1)'-Vv(:,1)/Sv(1,1)*Uv(:,1)'));
	%1/2*trace((S1-S2)(S2inv-S1inv))
	num2=.5*trace((Vv(:,1)/Sv(1,1)*Uv(:,1)'+Vc(:,1:1)/Sc(1:1,1:1)*Uc(:,1:j)')*(xval'-cal)*(xval'-cal)');
	%1/2*(S1inv+S2inv)(xso-xmean)'*(xso-xmean)
	DC(1,n+1)=abs(num1+num2);
end % if statement

end % function

function [srdin]= scalingvs(vs, scale,nMs)
% This function scales a cell array of matrices vs, across each row, by range scaling, rs,
% normalizing to max value, nmv, or by normalizing to unit legth,nul.
%   Note: Default scaing is nul unless otherwise selected.

if nargin<3, nMs=size(vs,1); end
if nargin<2, scale='nul'; end
srdin=cell(nMs,1);
for i=1:nMs
	srdin{i,1}=scalingvm(vs{i,1},scale);
end
end % function scaling

function [srdin]= scalingvm(vm, scale)
% This function scales a martix vm, across each row, by range scaling (rs)
% normalizing to max value (nmv) or by normalizing to unit legth (nul).
%   Note: Default scaing is nul unless otherwise selected.
if nargin < 2, scale='nul'; end
[n,m]=size(vm);
srdin=zeros(n,m);
if (strcmp(scale,'rs'))
	for i=1:n
		for j=1:m
			srdin(i,j)=(vm(i,j)-min(vm(i,:)))/(max(vm(i,:))-min(vm(i,:)));
		end
	end
else if (strcmp(scale,'nmv'))
		for i=1:n
			maxr=max(vm(i,:));
			for j=1:m
				srdin(i,j)=vm(i,j)/maxr;
			end
		end
	else if (strcmp(scale,'nul'))
			for i=1:n
				for j=1:m
					srdin(i,j)=vm(i,j)/norm(vm(i,:));
				end
			end
		else
			disp('Need to have a scalling method');
		end %(strcmp(scale,'nul'))
	end  %(strcmp(scale,'nmv'))
end %(strcmp(scale,'rs'))
end

function [ srdinred ] = inred(srdin, tuneR )
%This function  removes any rows, for the space merits in the in SRD input
%matrix, outside the tuneing parameter range. This is just for normal LOOCV


nMerits=size(srdin,1);
tuneRorg=tuneR;
if size(tuneR,1)==2
	tuneR=repmat(tuneRorg(1,:),nMerits-2,1);
	tuneR=[tuneR;repmat(tuneRorg(2,:), 2, 1)];
elseif size(tuneR,1)==1
	tuneR=repmat(tuneRorg, nMerits,1);
end

if nMerits==6 || nMerits==5 %just the X merits
	srdinred=cell(nMerits,1);
	srdinred{1}=srdin{1};
	for i=2:nMerits
		nRows=size(srdin{i,1},1);
		if isnan(tuneR(i-1,2)), 
			tuneR(i-1,2)=nRows; 
		end
		if (tuneR(i-1, 2)>nRows), 
			warning('The end tune is too big, the max possible is %g, for merit %g',nRows, i-1),
			tuneR(i-2,2)=nRows;
		end
		%assert(tuneR(i-1, 2)<=nRows,'The end tune is too big, the max possible is %g, for merit %g',nRows, i-1)
		assert(~(tuneR(i-1, 1)>nRows), 'The first tune is too big, the max possible is %g, for merit %g', nRows, i-1);
		nTuneNew=tuneR(i-1, 2)-tuneR(i-1,1)+1;
		for j=1:nTuneNew
			srdinred{i,1}(1:nTuneNew,:)=srdin{i,1}(tuneR(i-1,1):tuneR(i-1,2), :);
		end % for nTuneNew
	end % for nMerits
end % if nMerits==6


if nMerits==4 %just the X merits
	srdinred=cell(nMerits,1);
	srdinred{1}=srdin{1};
	for i=2:nMerits
		nRows=size(srdin{i,1},1);
		if isnan(tuneR(i-1,2)), tuneR(i-1,2)=nRows; end
		assert(tuneR(i-1, 2)<=nRows,'The end tune is too big, the max possible is %g, for merit %g',nRows, i-1)
		assert(~(tuneR(i-1, 1)>nRows), 'The first tune is too big, the max possible is %g, for merit %g', nRows, i-1);
		nTuneNew=tuneR(i-1, 2)-tuneR(i-1,1)+1;
		for j=1:nTuneNew
			srdinred{i,1}(1:nTuneNew,:)=srdin{i,1}(tuneR(i-1,1):tuneR(i-1,2), :);
		end % for nTuneNew
	end % for nMerits
end % if nMerits==4

end % function

