%Random matrix classification algorithm
%Alpha Lee 26/11/17
%
%Inputs
%
%training_binding: The training set for the binders 
%verification_binding: The test set for the binders 
%training_decoy: The training set for the decoys
%verification_decoy: The test set for the decoys 
%
%thres: a parameter that needs to be tuned such that the entire AUC curve
%is plotted (typically 100) 
%
%The datasets should be formatted as Nxp matrices, where N is the number of
%samples in the set and p is the number of descriptors per sample 

function AUC = RMD(training_binding,verification_binding,training_decoy,verification_decoy,thres) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%processing the binding training set 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tset_cleaned = training_binding; 

%compute z score 
[tset_cleaned_z, mu, sigma] = zscore(tset_cleaned); 
indzero_bind = find(sigma==0); %get rid of descriptors that have the same value for every member of the dataset  
tset_cleaned_z(:,indzero_bind) = []; 
mu(indzero_bind) = [];
sigma(indzero_bind) = []; 

%get covarience matrix and eigenvalues 
covar =tset_cleaned_z'*tset_cleaned_z/size(training_binding,1); 
[v, d] =eig(covar);

%Use the MP bound to get the number of significant eigenvalues  
p = size(tset_cleaned,2); 
n = size(tset_cleaned,1); 
l = diag(d);
num_eig = length(l(find(l>(1+sqrt(p/n))^2))); 

%get the significnt eigenvectors 
vv = v(:,end-num_eig+1:end); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%processing the decoy training set 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%get rid of columns with the same entry 
tset_d_cleaned = training_decoy; 

%compute z score 
[tset_d_cleaned_z, mu_d, sigma_d] = zscore(tset_d_cleaned);  

indzero_bind_d = find(sigma_d==0); 
tset_d_cleaned_z(:,indzero_bind_d) = []; 
mu_d(indzero_bind_d) = []; 
sigma_d(indzero_bind_d) = []; 

%get covarience matrix and eigenvalues 
covar_d =tset_d_cleaned_z'*tset_d_cleaned_z/size(training_decoy,1); 
[v_decoy, d_decoy] =eig(covar_d);

%Use the MP bound to get the number of significant eigenvalues  
p = size(tset_d_cleaned,2); 
n = size(tset_d_cleaned,1); 
l_decoy = diag(d_decoy);
num_eig_d = length(l_decoy(find(l_decoy>(1+sqrt(p/n))^2))); 

%get the significnt eigenvectors 
vv_decoy = v_decoy(:,end-num_eig_d+1:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%processing the binding verification set 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% first look at how close is the verification set to the binding set
verification_binding1 = verification_binding; 
verification_binding2 = verification_binding; 
verification_binding1(:,indzero_bind) = [];
verification_binding2(:,indzero_bind_d) = [];

%first look at how close the compounds are to the active training set 

%mean center and scale the verification set w.r.t. the active training set
veriset_mu = (verification_binding1-repmat(mu,size(verification_binding1,1),1))./repmat(sigma,size(verification_binding1,1),1);   
coeff = veriset_mu*vv; 

%project back into the ligand space 
proj_vect = (vv*coeff')';
norm_test = sqrt(sum((proj_vect-veriset_mu).^2,2));

%now look at how close the compounds are to the decoy training set 

%mean center and scale the verification set w.r.t. the decoy training
%set
veriset_mu = (verification_binding2-repmat(mu_d,size(verification_binding2,1),1))./repmat(sigma_d,size(verification_binding2,1),1);   
coeff = veriset_mu*vv_decoy; 

%project back into the ligand space 
proj_vect = (vv_decoy*coeff')';
norm_test_neg = sqrt(sum((proj_vect-veriset_mu).^2,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%processing the decoy verification set 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

verification_decoy1 = verification_decoy;
verification_decoy2 = verification_decoy;
verification_decoy1(:,indzero_bind) = [];
verification_decoy2(:,indzero_bind_d) = [];

%first look at how close the compounds are to the active training set 

%mean center and scale the verification set w.r.t. the active training set
veriset_d_mu = (verification_decoy1-repmat(mu,size(verification_decoy1,1),1))./repmat(sigma,size(verification_decoy1,1),1);   
coeff_d = veriset_d_mu*vv; 

%project back into the ligand space 
proj_vect_decoy = (vv*coeff_d')';
norm_test_decoy = sqrt(sum((proj_vect_decoy-veriset_d_mu).^2,2));

%now look at how close the compounds are to the decoy training set 

%mean center and scale the verification set w.r.t. the decoy training
%set
veriset_d_mu = (verification_decoy2-repmat(mu_d,size(verification_decoy2,1),1))./repmat(sigma_d,size(verification_decoy2,1),1);   
coeff_d = veriset_d_mu*vv_decoy; 

%project back into the ligand space 
proj_vect_decoy = (vv_decoy*coeff_d')';
norm_test_decoy_neg = sqrt(sum((proj_vect_decoy-veriset_d_mu).^2,2));

threshold = -thres:0.01:thres;  
for ii = 1:length(threshold) 
% compute false negative and false positve 

     true_pos(ii) = length(find(norm_test < (norm_test_neg + threshold(ii))))/length(norm_test);  
     false_pos(ii) = length(find(norm_test_decoy < (norm_test_decoy_neg + threshold(ii)) ))/ length(norm_test_decoy); 

end

AUC = trapz(false_pos,true_pos);

end